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The Physicist’s Companion to Current Fluctuations: 
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One of the main features of statistical systems out of equilibrium is the currents they exhibit in 
their stationary state: microscopic currents of probability between configurations, which translate 
into macroscopic currents of mass, charge, etc. Understanding the general behaviour of these cur¬ 
rents is an important step towards building a universal framework for non-equilibrium steady states 
akin to the Gibbs-Boltzmann distribution for equilibrium systems. In this review, we consider one¬ 
dimensional bulk-driven particle gases, and in particular the asymmetric simple exclusion process 
(ASEP) with open boundaries, which is one of the most popular models of one-dimensional trans¬ 
port. We focus, in particular, on the current of particles flowing through the system in its steady 
state, and on its fluctuations. We show how one can obtain the complete statistics of that current, 
through its large deviation function, by combining results from various methods: exact calculation 
of the cumulants of the current, using the integrability of the model ; direct diagonalisation of a 
biased process in the limits of very high or low current ; hydrodynamic description of the model 
in the continuous limit using the macroscopic fluctuation theory (MET). We give a pedagogical 
account of these techniques, starting with a quick introduction to the necessary mathematical tools, 
as well as a short overview of the existing works relating to the ASEP. We conclude by drawing the 
complete dynamical phase diagram of the current. We also remark on a few possible generalisations 
of these results. 
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I. INTRODUCTION 


Our world is a complex and chaotic place. In order to understand it better, we look for the fundamental 
laws that govern it, hoping that they are simple enough to be found, and universal enough to be useful. 
This is how, for instance, by observing the behaviour of various massive objects in different situations, 
Newton deduced his laws of motion, or how, by studying how gases react to changes in their environment, 
Clapeyron discovered the ideal gas law. Both of these laws are strikingly simple, and, although they are 
both approximations, they apply to an extremely wide range of systems. 

In fact, since a gas is a collection of massive objects, it must obey both laws, but at different levels of 
description: at the microscopic level, each individual atom follows Newton’s laws ; at the macroscopic 
level, the whole gas follows Clapeyron’s law. Moreover, Newton’s equation of motion, although concep¬ 
tually simple, has 6N variables for a gas with N atoms (the position and velocity of each atom), and 
the ideal gas law has only 3 (density, temperature, pressure). How these two descriptions are compatible 
is not entirely obvious, although this example is one of the simplest (but there are more complex ones, 
such as how neurons make a brain, to give an example from the other extreme). Understanding how one 
goes from one law to the other is just as essential as knowing the two laws themselves. It is the goal of 
statistical physics to bridge that gap, and to find out how simplicity can emerge from large numbers. 

Let us consider a system made of a large number of components, each obeying simple laws. There are 
several formalisms that we can use to that effect. If the system is isolated and stationary, with a well- 
defined interaction potential between its constituents, then it is reasonable to assume that its states are 
distributed according to the Gibbs-Boltzmann law: the system is at equilibrium. Given that assumption, 
one can then compute averages of macroscopic observables with respect to that distribution, and recover 
all the thermodynamic properties of the system, usually by obtaining the free energy, which gives the 
probability distribution of macroscopic variables. Non-analyticities in that free energy are the sign of 
phase transitions, which are the most interesting features of equilibrium systems, and have been studied 
extensively. However, in that framework, one is limited to static observables, as the Gibbs-Boltzmann 
distribution gives no information on dynamics. 

Should one be interested in dynamical observables, or if the system is driven by some external field and 
cannot be described by an interaction potential, one may take a step back and assume, not the form of the 
stationary distribution, but of the rates of transition between microstates and of the distribution of the 
time elapsed between successive transitions. This is equivalent to assume a certain probability density 
for dynamical trajectories in space and time, instead of a probability for a static configuration. The 
simplest and most common choice here is to assume Markovian dynamics: the evolution of a trajectory 
depend only on its state, and not on its history, which implies Poisson-distributed waiting times. Under 
reasonable assumptions, and if the dynamics do not depend explicitly on time, the system can be shown 
to relax to a steady state, the distribution of which is not a given, in contrast to equilibrium: it is a 
contraction of the distribution on trajectories, and is in general difficult to obtain. 

We may then distinguish two classes of Markovian dynamics: those that have detailed balance, and 
those that do not. In the first case, the steady state will have an equilibrium distribution, and the 
microscopic currents of probability between microstates will all be exactly zero: the state is not merely 
stationary, but entirely inert. Without detailed balance, some of those currents will be non-zero, and 
their magnitude is an indication of the work done by the environment in order to maintain that flow. 
That probability current, which translates into macroscopic currents of particles, charges, heat, etc. is 
what defines the system as being out of equilibrium. It is that type of systems which will be the focus 
of our attention. 

Non-equilibrium systems are quite ubiquitous in nature, and in particular in biology: any system 
which is meant to transport objects, cells, energy, etc. from one point to another is, by definition, out of 
equilibrium. In many cases, the system is one-dimensional (think of cells in a blood vessel, for instance, or 
molecular motors on actin), and the main quantity of interest is the flux that goes through it (see hgQ. 
The question is then to deduce the macroscopic behaviour of the system, and in particular the statistics 
of that current, from the microscopic definition of the model, and identify what, in that behaviour, may 
be generic among similar systems. 

There are many approaches to that problem. One of them, which has yielded major results in many 
other sub-fields of statistical physics, is to find a toy model which is simple enough to be mathematically 
tractable, yet complex enough to be physically relevant. In the case of equilibrium statistical physics, 
the Ising model has played such a role, and has had a central part in our understanding equilibrium 
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FIG. 1. Sketch of a one-dimensional non-equilibrium system. Particles move in a one-dimensional channel, 
between two reservoirs at fixed densities pa and pb- The particles interact with each-other, and are subject to 
a driving field V. Because of the field, and the possible imbalance between the reservoirs, there is a net flux of 
particles from one side of the system to the other. 


phase transitions. Its counterpart for non-equilibrium systems is the simple exclusion process: particles 
jump stochastically from site to site on a finite lattice, and may not jump to a site which is already 
occupied. One-dimensional versions of this model, with uniform jumping rates and either periodic or 
open boundary conditions, are exactly solvable, and allow for very precise calculations, such as relatively 
simple expressions for the full steady state distribution and current fluctuations. However, these results 
and calculations do not extend easily to other models, precisely because resolvability is not generic. 

Another approach is to start from an effective mesoscopic description, which requires more assumptions 
but applies to a broader range of systems. For equilibrium models, one may think, for instance, of the 
mean-field approach, where it is assumed that spatial correlations can be neglected. A similar framework 
for one-dimensional non-equilibrium particle systems is the macroscopic fluctuation theory, where the 
probability of a trajectory is assumed to depend only on the local averages of the density and current of 
particles. 

It is then quite natural to combine these two approaches: the first can be used to verify the assumptions 
made for the second, and the second can help in generalising the results obtained through the first. In 
this review, we propose to give an overview of the methods and results pertaining to these two approaches 
in the case of the one-dimensional asymmetric simple exclusion process with open boundaries. 

The layout of this review is as follows. 

In section [nj we give the main definitions and results related to the mathematical framework that is 
relevant to our interests: that of large deviations. In particular, we look at large deviations of time- 
additive observables for Markov processes in continuous time, and we define the so-called s-ensemble, 
which is a statistical ensemble for Markov processes where the current is seen as a free parameter. 

In section uni we get acquainted with the asymmetric simple exclusion process. In the first part of 
the section, we give the definition of the model and briefly review existing variants and results. In the 
second part, we look at the steady state of the system, first through the mean-field approach, and then 
through the exact expression for the stationary distribution (the so-called ‘matrix Ansatz’). 

In section EYi we present an exact expression for the complete generating function of the cumulants of 
the current in the open ASEP, and we take the limit of large sizes to extract the asymptotic behaviour 
of that expression, obtaining a different result for each phase in the ASEP’s diagram. We then look 
at what the corresponding behaviours are for the large deviation function of the current in each phase, 
which are valid for small fluctuations of the current. 

In section]^ we take the limit of infinitely high or low currents, and calculate the corresponding limits 
for the large deviation function through direct diagonalisation of the conditioned process. In the low 
current limit, we obtain a perturbative expansion around a diagonal matrix. In the high current limit, 
we obtain a system equivalent to an open XX spin chain, diagonalisable through free fermion techniques. 

In section |Vlj we use the macroscopic fluctuation theory to obtain the full dynamical phase diagram of 
the current for the open ASEP. We compare these results with those obtained from the exact calculations 
of the previous sections. We also remark that the method used in that section is in principle applicable 
to any one-dimensional bulk-driven particle gas. 
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The present review is largely based on the author’s PhD manuscript [I]. It is intended to give a 
self-contained and reasonably detailed account of the tools involved in determining the large deviations 
of the current in the asymmetric simple exclusion process, as much as of the results that they yield. 
For the sake of brevity and legibility, some of the finer details of those calculations, as well as most 
technical aspects of everything related to the integrability of the model, have been omitted, but should 
the readers be in need of clarifications, they may refer to [T] or to the references given in section 
That section contains most of the bibliographical references of this review, and although it is far from 
being exhaustive, it should provide an adequate starting point for the curious reader. 
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II. A CRASH COURSE IN LARGE DEVIATIONS 


This first section contains a brief introduction to the mathematical objects that we will be manipulating 
in the rest of the review, namely: large deviation functions. We first define the concept of large deviations 
in general, and give a few useful theorems. We then apply the concept to time-additive observables for 
Markov processes in continuous time, and look at two specific examples with interesting properties: the 
time-integrated empirical vector, and the entropy production. For a thorough review of this topic, one 
may refer to [5] and [3]. 


A. Definition and a few useful results 


Consider a system defined by a size N, and an observable a intensive in N, which has a probability 
distribution P;v(fl) for each N. It is said that a obeys a large deviations principle with rate g{a) if the 
limit 


ff(a) 


lim 

N—>-oo 


log(Pjv(a))- 

N 


( 1 ) 


is defined and finite for every a. In other terms, g{a) is the rate of exponential decay of PAr(a) with 
respect to N. Its minimum is the most probable value of a. We then write: 


PAr(a) « 


( 2 ) 


where the « signifies precisely what is written in eq.Q. 

Note that N is not necessarily an actual size or a number of elements: it can be a time span, a number 
of events, or any variable that can be taken to infinity. Also note that a doesn’t have to be a scalar 
observable. It can be a function (in which case g[a] is a large deviation functional), or any mathematical 
object for which a probability can be defined in the system under consideration. 


From the large deviation function of a, one can obtain that of an observable b = f{a), where the 
function / is not necessarily injective, by writing 

f y V f TiT /■ \ —A/” min [g(a)l 

P 7 v(&) = y PAr(a)d(&-/(a))da = y /(a))da ~ e , (3) 


the last expression being obtained from a saddle-point approximation for large N. If / is injective, one 
simply obtains a change of variables from a to b. If, on the contrary, / is not injective, which is to say 
that 6 is a contraction of a, we get the contraction principle: the large deviation function g(b) is given 

by 


Kb) = min [g(a)]. 

f{a)^b 


(4) 


In principle, a may be any mathematical object and / any function, but we will mostly consider linear 
transformations, where a is a vector and / a matrix with a non-maximal rank. 


An alternative way to treat a probability distribution which decays fast enough is through its rescaled 

OO 

cumulants Fife, or their exponential generating function E{p) = defined through 

k=l 


= J P^(a)e'^^“da (5) 

which is to say that the generating function of cumulants is the logarithm of the exponential generating 
function of moments. 

If we replace PAr(a) by its limit under the large deviation principle: PAr(a) —>■ we get, in the 

large N limit: 

^NE(p) f ^-N{g(a)-t,a)^^^ (-g) 
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which yields, through a saddle-point approximation, 


E{^) = max[/ia — 5 (a)] 

a 


(7) 


or equivalently 


E{fi) = fia* - g{a*) 


da 


g{a*) = g 


( 8 ) 


where a* is the value of a at which the maximum in eq.Q is attained. That is to say that E and g 
are Legendre transforms of one another and contain essentially the same information (unless g has a 
non-convex part). The inverse transformation formula is then 


5 ( 0 ) = min[/ra — E{p)] 


(9) 


or 


g{a) = g*a - E{g*) , —E{g*) = a 


( 10 ) 


where g* is the value of g at which the maximum in eq.( 10 1 is attained. This last equation is part of the 
Gartner-Ellis theorem, which states that if E{g), defined through eq.®, is well-behaved, then a obeys 
a large deviation principle with a rate g{a) obtained through eq. § PT 


One of the most useful features of the cumulant approach is how it combines with the contraction 
principle. Consider a vector a and a non-injective matrix /. The generating function of the cumulants 
of a contracted observable b = f ■ a is given by 

E{g) = max[/i • b — g{b)] = max[/i • b — min [ 5 (a)]] = max[/l ■ f ■ a — 5 (a)] = E{jl ■ /) (11) 

b b fa=b a 

which is to say that the function E is in fact the function E applied to a variable g ■ f which has fewer 
degrees of freedom than g, (because g is conjugate to b which has fewer degrees of freedom than a). In 
other words, contracting at the level of cumulants reduces to taking special values of g, which is often 
much easier than finding the minimum in eq. @- 

One last remark we should make is that P 7 v(a) bas a sub-exponential pre-factor which we may (and 
did) neglect entirely, as long as it has no poles in a. In the case that it does, all the saddle-point 
approximations that we performed have to be modified to take into account contour integrals around 
these poles, which may dominate the large N limit [1]. 


B. Dynamical large deviations for Markov processes 

We will now see what can be said of large deviations in the context of continuous time Markov processes 
on a finite state space. 


1. Definitions 

Consider a Markov matrix M acting on states {C}, with rates w{C',C) from C to C, and escape rates 

^(C) =J2w{,C',C): 

C 


M=J2 wiC,C')\C){C'\ -J2riC)\C){C\. 

C/C' c 

This defines the time evolution of a probability vector \Pt) through the master equation 


( 12 ) 




( 13 ) 
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of which the solution is, formally: 


|Pt)=e‘"|Po) (14) 

where |Po) is the initial probability distribution. In the limit of long times, if M is not reducible, this 
vector converges to a steady state 

lim \Pt) = \P*) with M\P*) = 0. (15) 

t—¥oo 

Equivalently, we can define the probability density P(C(t)) of a history C{t) going through configu¬ 
rations Ci with waiting times tf. the Markovianity of the process tells us that the waiting times are 
Poisson-distributed with respect to the escape rates, which gives us 

P[C(t)] = e-‘~’’(^~)rt;(CAr,Civ-i) (16) 

where time increases from right to left. This can then be recast in a more compact form: 


N-l 

log(p[C(t)]) = - r{C{t))dt+ ^ logw;(C*+i,C*). 

do 


(17) 


The entries of e‘^^ are then given by the sum of the probabilities of all paths that start and end at the 
corresponding microstates. 


2. Time-additive observables 

We can now define time-dependent observables on those histories and look at their large deviations in 
the long time limit. Let us consider an observable At defined as a functional of a history C{t): 

^=F[C(t)]. (18) 

For the sake of simplicity, we are only interested in observables that are additive in time, meaning 
that if a history C{t) is the concatenation of two shorter ones Ci{t) and C 2 (t), which we will write as 
Ci{t) 0 C 2 {t), the functional F distributes over them: 

F[Ci{t)(BC2{t)]=F[Ci{t)]+F[C2{t)]. (19) 

This forces F to be local in time (independent of time correlations). We also assume F to be time- 
invariant, i.e. independent on the value of the initial time to in C(t). 

These constraints allow us to find the general form of such functionals. Consider first a history without 
any transitions: C(t) = Ci. In this case, time additivity can be used to show that F[C{t)] is proportional 
to the duration ti of the process. The proportionality coefficient may depend on Ci, and we will call it 
V (Cl). Consider now a history with one transition: the system is in Ci for a duration ti, and in C 2 for a 
duration f 2 - By cutting out, using additivity, the portion of history before ti — e and that after ti 0 e, 
with e going to 0, one is left with just the transition, which has a contribution to F that depends only 
on Cl and C 2 . We will call it [/(C 2 ,Ci). 

Putting those pieces together, and considering that any history can be decomposed into portions 
containing at most one transition, we can finally write: 


F[C(t)] 



N 

V{Cit))dt + Y,U(C,,C,.i) 

Z=1 


( 20 ) 


which is expressed schematically on fig.-[^ 

The first part of that expression, containing P, is a state observable, depending only on the empirical 
vector, which is to say the relative time spent in each microstate C. It contains no direct information 
on the currents between microstates, although it may depend indirectly on the system being in or out 
of equilibrium, through time-correlations. The second part, containing U, is a jump observable, and is a 
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FIG. 2. Functional F over a schematised history. Each straight portion contributes a simple term to the whole 
function. Waiting periods (in blue) give a contribution that is extensive in time, and depends only on one 
configuration. Transitions (in red) give a term that depends on the two configurations involved. 


direct measure of those currents. Each U{Ci,Ci-i) can be seen as a counter for the transition between 
Ci-i and Cf. every time it is used in the evolution of the system, the value of At increases by one quantum 
of U. Whether each value of U is taken as an independent variable, or given a precise value, determines 
what information is monitored regarding the way those transitions are used, but in many cases, one thing 
that makes a crucial difference on the resulting behaviour of At is whether the system is in equilibrium 
or not, as we will see shortly. 


But first, let us have a look at the cumulants of this observable. The generating function Et{v) of the 
cumulants of at = \At (the intensive version of At) can be expressed as: 




(e‘"“‘) =y’ e"’'f^(*)lp[C(t)]I?[C(t)] 


( 21 ) 


where 'D[C{t)\ is the measure associated with histories (this can be simply defined in the discrete time 
case, and then taken as a formal limit for small 5t). 

Replacing P[C(t)] by the expression in eq.([T7|), we see that 


log(e"F[C(‘)lp[C(t)]) 


j, N-l 

J (r(C(<))-7^R[C(t)])dt+^ (logu;(C,+i,C,) + 7^C/(C,+i,C,)) 




( 22 ) 


which is to say that it is the (un-normalised) weight produced by a modified Markov matrix Mj^ defined 
as 


- J2(^iC) - EViC))\C){C\ 

c^c c 


(where we take U{C,C) = 0). By replacing M by in ( [l4| ), we get 

|P.(t)) = e*^lPo) = E / |C^) 

Cn 


(23) 


(24) 


where, as intended, the probabilities of histories have received an extra factor. We can finally 

sum over the final configuration Cjv by projecting to the left on the uniform vector (1| (of which all 
entries are 1) and write 


= (l|P,(t)) = (l|e*^-'|Po). 


(25) 
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Note that we use u as a, generic parameter, which can in fact be a function of the configuration (as in 
section 


IIB 31 or of the transition. 


Moreover, as long as v, V and U are real, the Perron-Frobenius theorem applies: the largest eigenvalue 
Ai, of is non-degenerate. If we write the corresponding eigenvectors as \Pu) and (Pi,|, we get, for 
large times, 






(26) 


By combining equations (25) and (26) for t large, we get « e‘^''(I|Py)(Pi,|Po), which is to say that 


E{v) = A^. 


(27) 


The generating function of the cumulants of any additive observable is therefore equal to the largest 
eigenvalue of the associated deformed Markov matrix M^. This is a classic result from the Donsker- 
Varadhan theory of temporal large deviations [5H9]. Notice that this is a property of the deformed 
matrix, and not of the initial or final configurations: regardless of those, the long time behaviour of the 
generating function of the cumulants will be the same. One can easily make sense of this: if the duration 
of the process is large enough, then the system only takes a small time at first to reach its steady state, 
and gets out of it very near the end. The rest of the evolution can be considered to be around the steady 
state, whatever the initial and final distributions are, and the part of At which is extensive in time comes 
only from there. The initial distribution only gives a time-independent term (l|P,y)(P,^|Po), which is 
negligible (unless, as we mentioned earlier, it has poles). 


The eigenvectors \P„) and {P^\ also carry information on current fluctuations. Considering eq.(24) for 
t large, we can regroup all the histories for which F[C(t)] = tf and the final configuration is Ct- Writing 
that probability P(/ & Ct) as P(Ct|/)P(/) (where P{A\B) is the probability of A, conditioned on B), 
and invoking the large deviations principle P(/) ~ we get 


e‘^M|P,)«^|Ct) f P(Ct|/)e‘("/-s(/))d/. 


(28) 


Finally, just as in eq.([^, a saddle-point approximation on yields and fixes the value of 

/ to ^P(z/). Injecting this in (28), we get 


PAC)=P{Ct=C I f=-^Eiu)). 


(29) 


This tells us that the vector \P^) is in fact the probability vector of the final configuration, knowing that 
the value of / through the evolution of the system was ■^E{v) [10]. A similar calculation on the left 
eigenvector (P,,| shows it to be the probability vector of what the initial configuration was, knowing that 
the value of / is ^Eiv): 


Pu 




(30) 


Finally, the product of the two gives the probability of observing a configuration at any point during the 
evolution of the system (but far enough from the initial and final time), conditioned on the value of /: 



d \ 

P,{C)P,iC) = p{c 



(31) 


Note that these vectors are implicitly normalised appropriately. 

Those three distributions are quite different from one another, and in particular their most probable 
states are in general not the same. We will be mostly interested in eq.(311, for reasons that will become 
apparent later, in section IV When the observable of interest is the entropy production (which we will 
examine in section IIB 4), the statistical ensemble defined by those probabilities, where v is considered as 
a free parameter (similar to a temperature, as it is conjugate to the entropy), is sometimes called the ‘s- 
ensemble’ m- It is quite natural to consider this ensemble: since entropy production plays an important 
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part in the system being out of equilibrium, being able to control its value and look at how the system 
responds provides us with useful information on its behaviour. Moreover, much like in equilibrium, 
where one can think of a Lagrange multiplier as an extra interaction, one may think of as equivalent 
to a process with modified dynamics, which has typical values for observables given by the fluctuation 
values in the original process nmg. The Markov matrix for that new process can be obtained through 
a so-called ‘Doob transform’, which consists in shifting its eigenvalues by —E{v) and conjugating the 
matrix through so that its largest eigenvalue be 0 and its dominant left eigenvector be uniform. Its 
steady state is then given by Py(C)Py(C), which is another argument for that distribution being the most 
natural one. 

In the following sub-sections, we will consider two special cases of time-additive observables: the 
empirical vector, which is a state observable, and the entropy production, which is a jump observable. 


3. Large deviations of the empirical vector 


The first special case that we will consider is that of the time-integrated empirical vector, which is the 
vector of the total time spent in each configuration. We will write its intensive counterpart as p, defined 

by 


1 /■* 

P(c) = 7 / Scfi(,T)dT 

b Jo 


(32) 


which contains the fractions of time spent in each configuration c. We recognise a special case of eq. (20) 


where V{C) = 6c,c and U{C',C) = 0. This is the most general state observable, which can give any other 
through contraction. 

The deformed Markov matrix corresponding to p is given by 


w{C,C')\C){C'\ - ^(r(C) -hc)\C){C\ = M + H 
C/C' c 


(33) 


where P is a diagonal matrix with entries he- The largest eigenvalue E{h) of contains all the 
cumulants of the empirical vector, including n-point functions. For instance, 

m = {piCi)p{C2))t - {p{C0)t{piC2))t (34) 

d/icid/ic2 

where (•)( refers to the average in time (which we will soon need to distinguish from an ensemble average). 


There is another empirical vector tt which we may define, by looking at the probability of being at C 
at a given time t, averaged over a large number N of copies of the process. If the probability vector at 
time t is |Pt), the exponential generating function of the cumulants of tt is given by 

£{h)=log{{l\e^\Pt)) (35) 

which involves ensemble averages {-je- 

It is important to note that these two generating functions, as well as the corresponding large deviation 
functions, are in principle different: fluctuations over time are not identical to fluctuations among copies 
of the system. They do however have one common feature. We may write the first cumulant from E{h) 
as 




= lim - 
h=0 t 


P^(C)dT = P*(C) 


(36) 


where we recall that P* is the distribution of the steady state of M. That same first cumulant from S{h) 
gives, in the limit of long times: 






h=0 


= P*(C). (37) 


Combining this identity with any contraction of p or tt, we conclude that time averages and ensemble 
averages of state observables are identical, which is to say that the Markov process we are considering 
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is ergodic. All the other cumulants are in principle different, which can be understood by the fact that 
time-correlations (i.e. the dynamics of the system) play a role in E which they don’t in £. This implies 
that the absolute minima of both large deviation functions have the same locus, but that their shape 
(stiffness, skewness, etc.) around those minima is different. 

It is also interesting to note that this remains true for averages in a conditioned process. Let us for 
instance take a generic deformed Markov matrix with dominant eigenvectors |P^) and (P^|, and 
consider the two generating functions 

Ei^i, h) = lim I log((l|e‘('^-+") |Po)) (38) 


and 


£{fi,h)=log{{P^\e^\P^)) (39) 

(where the product between the two eigenvectors \P^) and (P^| means we are looking at a time which is 
somewhere in the middle of a long evolution). A calculation similar to what we did earlier yields: 




d 

dhc 


A(m, 0) = P^(C)P^(C) 






(40) 


where the product P^(C)P^(C) is normalised. 


4- Large deviations of the entropy production 

We now consider a special jump observable: the entropy production St for an evolution between times 
0 and t, defined as 


S.|C(r)|=log(jpi), (41) 

which measures how probable a history C(r) is compared to its time-reversal C^{t) = C{t — r). It 
corresponds to a time-additive observable with 

P(C)=0 , t/(C,C')=log(u;(C',C))-log(u;(C,C')). (42) 

so that the corresponding deformed Markov matrix is given by: 

= ^w{C,C'f+''w{C' ,C)-''\C){C'\. (43) 

c.c 

We notice that it has a peculiar symmetry: if we replace v hy —l — v, the exponents in are swapped, 
which has the same effect as transposing the matrix. That is to say. 




(44) 


This has interesting consequences on its eigensystem: all the eigenvalues are symmetric with respect to 

v —l — v, and the associated right and left eigenvectors are exchanged. This is the famous ‘Gallavotti- 

Cohen symmetry ’ [HHm, which is one of the only universal features known for non-equilibrium systems. 

In particular, the generating function of the cumulants of the intensive entropy production s, and the 
conditional probabilities that we have defined earlier, verify: 

E{iy)=E{-l-^) , P,{C) = (45) 

which becomes, for the large deviation function, 

9{s) - gi-s) =-s (46) 
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or, equivalently, 


P(-s) =e-‘"P(s) (47) 

for t —>■ oo. This last equation is called the ‘fluctuation theorem’. It was first observed by Evans, Cohen 
and Morriss in m, then proven by Evans and Searles in m, and later led to Gallavotti and Cohen’s 
formulation of their symmetry. The theorem means that a negative entropy production rate is much less 
probable than its positive counterpart, but not impossible. This does not, as it might seem, contradict 
the second law of thermodynamics, which is expressed only for the average of s. It in fact validates it, 
since it implies that the mean value of s must be positive. 

We may finally note that, for equilibrium systems, where the detailed balance condition imposes that 

P*{C)w{C',C) = P*{C')w{C,C') (48) 

for any two configurations C and C', the deformed Markov matrix turns out to be similar to the un¬ 
deformed one, so that they have the same eigenvalues. Consequently, the generating function of the 
cumulants of the entropy production rate is identically zero, and its large deviation function is a delta 
function: 


E{v) = 0 and P(s) = ^(s). 


(49) 


There is therefore no entropy production whatsoever in the case of an equilibrium system. This, as in 
eq.(27), is a property of the transition rates, and not of the initial or final distributions. There can still 
be a conservative exchange of entropy between the initial and final configurations of a history (if their 
equilibrium probabilities are not equal). 
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III. THE ASYMMETRIC SIMPLE EXCLUSION PROCESS 


In this section, we are introduced to the Asymmetric Simple Exclusion Process. After giving its 
dehnition, we briefly go through the existing literature related to that model, including its many variants 
and connections with other problems in physics, mathematics and biology. We then examine the steady 
state of the model, first through a mean-field approach and then through an exact calculation. 


A. Definition of the model and variants 

Consider a one-dimensional lattice with L sites (or a row of L boxes), numbered from 1 to L. Each 
site can be empty, or carry one particle. Those particles jump stochastically from site to site, with a rate 
p if the jump is to the right, from site i to site z -I-1 (and which will be set to p = 1 by choosing the rate 
of forward jumps as a time scale), and a rate g < I if the jump is to the left, from site i to site i — 1. 
The jumping rate is larger to the right than to the left in order to mimic the action of a field driving 
the particles in the bulk of the system. Each end of the system is connected to a reservoir of particles, 
so that they may enter the system at site I with rate a or at site h with rate i5, and leave it from site 1 
with rate 7 or from site L with rate /3. Those rates allow us to define the effective densities of the two 
reservoirs. In all of these operations, the only constraint that must be obeyed is that of exclusion, which 
is to say that there cannot be more than one particle on a given site at a given time, so that a particle 
cannot jump to a site that is already occupied. These rules are represented schematically on hg.-[^ 


OL 1 

o, n 



u uvv ^ 

7 ^ 5 


FIG. 3. Dynamical rules for the ASEP with open boundaries. The rate of forward jumps has been normalised 
to 1. Backward jumps occur with rate q < 1. All other parameters are arbitrary. The jumps shown in green are 
allowed by the exclusion constraint. Those shown in red and crossed out are forbidden. 


Configurations of the system are written as strings of O’s and I’s, where 0 indicates an empty site, and 
1 an occupied site. The Markov matrix governing that process is a sum of local jump operators, each 
carrying the rates of jumps over one of the bonds in the system: 


with 


L-l 

M = Too + ^ Mi -I- 

i=l 





■0 

0 

0 

O' 




—a 7 

, M, = 

0 

-q 

1 

0 


'-6 13' 

mo = 

a —7 

0 

q 

-1 

0 

, rriL = 

s 




0 

0 

0 

0 . 




(50) 


(51) 


It is implied here that mg acts as written on site 1 (and is represented in basis {0,1} for the occupancy 
of the first site), and as the identity on all the other sites. Likewise, ttil acts as written on site L, and 
Mi on sites i and i + 1 (and is represented in basis { 00 , 01 , 10 , 11 } for the occupancy of those two sites). 
Each of the non-diagonal entries represents a transition between two configurations that are one particle 
jump away from each other. 

As we mentioned earlier, because of the asymmetry of the jumps, the particles flow to the right, which 
results in a macroscopic current deeply connected to the non-equilibrium nature of the system. It is a 
special case of time-additive observable (as presented in section IIB 2) where V = 0 and U is taken as 
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1 for the transitions where a particle jumps to the right, and —1 for those where one jumps to the left. 
We will see in section IV that it is strongly related to the entropy production (they are in fact equal, up 
to a constant), and we will be referring to probabilities conditioned on the current as the s-ensemble as 
well. All the results presented in this review are related to the characterisation of that current and its 
fluctuations. 


1. Variants of the A SEP 

There are a few simpler cases that one can consider. The first is to force the particles to jump only 
to the right, by taking q = j = 6 — 0. In this case, the model is called the totally asymmetric simple 
exclusion process (or TASEP), and we will often use it in our calculations, as its behaviour is identical 
to that of the ASEP for all intents and purposes, but much easier to deal with. 

The second is the opposite limit, where the jumps are as probable to the left as they are to the right: 
<7=1. This case is called the symmetric simple exclusion process (or SSEP). It is an example of a 
‘boundary-driven diffusive systems’ (as opposed to the ASEP, which is bulk-driven by the asymmetry, 
and is therefore not diffusive). Its behaviour is quite different from that of the ASEP, and we will not 
consider this limit in the present review. 

Sitting somewhere between the SSEP and the ASEP is the weakly asymmetric simple exclusion process 
(or WASEP), where the asymmetry 1 — g is taken to scale with the size of the system as L~^. This 
is done in order to make the integral of the field in the bulk, which is of order L(1 — q), comparable 
with the difference of chemical potential between the reservoirs, which is a constant with respect to L. 
The ASEP and the WASEP correspond to two different ways to take the large L limit in the system: in 
the ASEP, no rescaling is done to the driving held, so that the large size limit corresponds to a system 
of increasing length, with the lattice spacing remaining constant, which is relevant to model a system 
which is really discrete (think for instance of ribosomes on a long string of mRNA, or any other example 
of discrete biological transport). In the WASEP, on the contrary, the held is rescaled as L~^, so that 
the large size limit corresponds to a system of hxed length, with a smaller and smaller lattice spacing, 
going to a continuous system when L reaches inhnity. We will be using the WASEP as a starting point 
in section ED 

One can also consider different geometries for the model. Take for instance the ASEP with periodic 
boundary conditions, i.e. on a ring (hg.-|^b). In this case, the system is not connected to any reservoir, 
and the number of particles is conserved. This makes it somewhat easier to deal with: the steady-state 
distribution is uniform, and the coordinate version of the Bethe Ansatz can be used to solve it, as we 
will see in section HVB II 

The ASEP can be dehned on an inhnite lattice instead (c.f. lower part of hg.-|4 d). In this case, there 
is in general no convergence to a steady state (for generic initial conditions), and the observable of choice 
is instead the large time behaviour of the transient regime. 

Finally, one can put more than one type of particles in the system, and consider the multispecies ASEP 
(%•§ c). The exchange rates must then be defined between any two different species of particles. The 
simplest case to consider (and the most tractable one) is that where the types of particles are numbered, 
from 0 (for holes) to K (for the ‘fastest’ particles), and where a particle of type k sees all lower types 
k' < k as holes, which is to say that the rates of exchange of two particles of types ki < k 2 are 1 for 
^ 2^1 — kik 2 and q for kik 2 k 2 ki (those rates are represented on fig.-|^c, where different species of 
particles bear different colours, and are numbered by their rank). 


2. Brief overview of the ASEP’s family tree 

We mentioned biological transport earlier for a good reason: the first definition of an ASEP-like model 
was made in 1968 in uniiin] precisely in order to study the dynamics of ribosomes on mRNA (fig.-|^a) . It 
is still used today in that context, often with a few modihcations to make it slightly more realistic, such 
as making the particle reservoirs finite m or even shared between several systems [52], changing the 
jumping rates from site to site |23|, changing the jumping cycle by adding an inactive state for particles 
[24| . allowing them to attach or detach in the middle of the chain jSS], and so on. These are only a few 
recent examples, but a thorough review can be found in |26j . 
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It has also been noticed that the ASEP is strongly related to the XXZ spin chain with spin ^ m- 
the Markov matrix of the ASEP and the Hamiltonian of the spin chain are related through a matrix 
similarity. This fact goes deeper than a simple mapping between two systems: the XXZ spin chain is well 
known and well studied, as it has the mathematical property of being ‘integrable’, meaning that it can 
be solved exactly, for instance through the Bethe Ansatz |5H], and that we can expect precise analytical 
results from it [29]. For that reason, many results have been obtained for the ASEP by adapting the 
Bethe Ansatz to its formalism |5DH55] . and even results that have been found by other means (such as 
those presented in [33]) are in fact consequences of that property. The downside of this fact, one might 
argue, is that those methods are only transposable to other integrable systems, but the undeniable upside 
is that we might have access to very precise results, which could lead to discovering universal features of 
non-equilibrium systems. 


Another model related to the ASEP is that of random surface growth [ID] (fig-j^d). In that model, 
a wall, made of square blocks (with corners pointing up and down), grows by a procedure where blocks 
can fall in valleys at rate 1, or lift off of peaks at rate q. The relation to the ASEP is rather obvious if 
one replaces upward slopes (when reading from left to right) by holes, and downward slopes by particles 
(adding a block means replacing ‘down up’ by ‘up down’, i.e. 10 by 01, and removing one is the opposite 
operation). In this context, the situation that is usually considered is that of the infinite ASEP, with a 
simple initial condition, such as a given mean density to the right of site 0 and another one to the left 
(the simplest one being all Is to the left and all Os to the right [H] 32] : as represented by the dashed line 
in fig.-|^d), although more general ones can be considered |43j . One of the most interesting quantities, 
just as for finite size, is the total current of particles that went over the bond at the centre of the system, 
which is equal to the number of blocks that have been added above that site, i.e. to the height of the 
surface (represented in green in fig.-|^d; each green block corresponds to one of the particles that crossed 
to the right of the system). After a first breakthrough by Johansson in [IT], the fluctuations of that 
height were conjectured [44] and then proven [45] to be related to the famous Tracy-Widom distributions 
governing the eigenvalues of random matrices [lillT]. Some even more complex quantities have been 
studied, such as the general n-point correlations of the height [35]. Moreover, there is a whole class 
of systems, the KPZ universality class (named after Kardar, Paris! and Zhang, authors of the seminal 
article that started it all [40]), that are governed by the same laws, such as the directed random polymer 
gaisn], or the delta-Bose gas m- One can in particular recover the KPZ equation [S2] from the WASEP, 
or an equivalent solid on solid model, through a particular rescaling [S3], as well as some of the critical 
exponents through renormalisation group techniques [54]. 

One of the main features of the KPZ universality class is its dynamical exponent 2 = |: in the 
language of the exclusion process, the time it will take on average for a local perturbation of the density 
of particles, or a marker particle, to spread over a region of size L scales as t ^ as opposed to 

t ^ for a diffusion and t L for a ballistic system. That property relates to the relaxation of the 
system at long but finite times, which we will not consider in the present work (as we will be looking 
at steady states, i.e. at infinite times), but we will observe the appearance of that exponent in certain 
places, such as in equation (166). 

Experimental evidence of the relevance of that model has been obtained in a liquid crystal undergoing 
a phase transition [55]. This subject has generated many more works than could be summed up here, 
and the reader can find more information in reviews such as 


The ASEP can be related to many more models and mathematical objects, such as chains of quantum 
dots [60], alternating sign matrices [6ll[^ (through its connection to the XXZ chain), continued fractions 
[63] , Brownian excursions [64U66] , Askey-Wilson polynomials [671 - 169] , and a large family of combinatorial 
objects which all have a connection to Catalan numbers IZD]. 


3. Earlier results 

All these interesting connections notwithstanding, the ASEP is a very popular model in itself [IIHZS] 
(it has even been referred to as the Ising model of non-equilibrium systems [74]), and has been the subject 
of a tremendous number of works. 

The SSEP, for one, has established itself as an archetype of diffusive systems with interactions, for 
which many universal results have been found, such as the cumulants of the current in a periodic system 
m or an open one m- Those results all have to do with the so-called ‘macroscopic fluctuation theory’ 
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FIG. 4. The ASEP’s family album: a) Ribosomes on mRNA, b) Periodic ASEP, c) Multispecies ASEP, d) Surface 
growth. 
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(or MFT) [TTHSO] . developed to deal with the fluctuations of diffusive systems through a hydrodynamic 
approach [ST]. As for results more specihc to the SSEP or the WASEP, the large deviation functional 
of the density profiles was expressed in [55], leading to the joint large deviations functional for the 
current and the density |83| which we will be using in section VI The cumulants of the current for the 
open SSEP were found in |84] . and were observed to depend on a single variable and not on the two 
boundary densities independently. This lead to the discovery of a surprising symmetry connecting the 
non-equilibrium SSEP (with different reservoir densities) to a system at equilibrium [55H57] . The full 
cumulants of the current for the periodic WASEP were found in [55]. In that case, as was found in 
[55] and further analysed in EH! El], the system undergoes a phase transition in the s-ensemble where, 
for a low enough current, the optimal density profiles become time-dependent. A similar transition 
can be found for the activity (the sum of jumps, regardless of direction) in the SSEP |92j . or even for 
combinations of current and activity |93j . Recently, the large deviations of the position of a tracer in the 
one-dimensional single file diffusion (a continuous equivalent to the SSEP) using MFT [SUES]- See [75] 
for a review of some of these results. 


The periodic ASEP, with its hxed number of particles and its trivial steady state (all the configurations 
are equally probable, as long as they have the correct number of particles), has mostly been studied for 
the fluctuations of the current. The full generating function of those was found for the TASEP in 
1998 EHIET], and although the second cumulant for the ASEP was found prior to that [55], the complete 
generating function was only obtained more than 10 years later [5DH55]. Some other results were obtained 
for the periodic TASEP, such as the gap (i.e. the characteristic time of the transient regime) [551 llOOj . 
and, very recently, the whole distribution of the spectrum of the Markov matrix and the behaviour of 
the low excitations [101H103] . The s-ensemble was also investigated, for the limit of very large currents, 
and the probabilities of the configurations were found to be those of a Dyson-Gaudin gas (the discrete 


analogue of a Coulomb gas) [104] . We will come back to that last observation in section V B 


The open ASEP is richer than the periodic case, but much harder to handle. The structure of the steady 
state itself is quite intricate: it was first found in m for the TASEP thanks to some surprising recurrence 
relations between the weights of the configurations for successive sizes. It was then generalised to the 
ASEP by expressing those relations in algebraic form [55], giving birth to the ‘matrix Ansatz’, which we 
will present in section IIIB 2[ Depending on the values of the two reservoir densities, the system can find 
itself in three different phases, which was discovered for the TASEP in [106] as an interesting feature of 
non-equilibrium systems (since, for equilibrium systems with short range interactions, transitions cannot 
be induced by boundaries). This phase diagram was refined in |107j where sub-phases were found with 
different correlation lengths. Those results were extended to the ASEP in [53 1108] (part of which we 
present in section III B 2\ . The 2-point correlation function |109j and then the complete n-point function 
[55] were calculated for the TASEP, for some values of the boundary densities, and the same was later 
done for the ASEP in mm- Most of these results rely on the matrix Ansatz, and a review of those 
results and methods can be found in mu. See also m for a review of various results for the steady 
state of the open ASEP. More recently, a matrix Ansatz was found for the steady state of the two-species 
open ASEP [I12j . 

Other properties of the steady state were analysed, such as the static density, current and activity 
distributions [11311114] , the large deviation function of the density profiles [115] , or the reverse bias regime 
(where the boundaries impose a current floving to the left) [11611117] . A hydrodynamic description, 
named ‘domain wall theory’ (or DWT) [ST] 111811121| where states of the system are approximated by 
regions of constant density separated by discontinuities called shocks, was proposed to describe the large 
scale dynamics of the system, even in the transient regime, and the hydrodynamic quasi-potential was 
obtained along with the optimal relaxation pathways in but no full equivalent of the MFT has yet 

been devised. 

One of the main reasons why the open ASEP is more difficult to study than its circular sibling is that 
the Bethe Ansatz cannot be used as easily in this case. The coordinate version of the Ansatz which 
we presented in section [III B 1[ where the particles are treated as plane waves, relies on the number 
of particles being fixed, and breaks down in the open case. Variants of the coordinate Bethe Ansatz 
have been used successfully to build excited eigenstates of the system for some special cases of the 
boundary parameters [3411351 155] (in particular with triangular boundary matrices), and, in conjunction 
with numerical analysis, to find the relaxation speed of the system (i.e. the gap of the Markov matrix) 
[5511571 [T55]. as well as the asymptotic large deviation function of the current inside of the Gaussian 
phases [124] . A generalisation of the matrix Ansatz was used in [125] to calculate the second cumulant of 
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the current in the totally asymmetric case. In |126j . an alternative to the Bethe Ansatz, named the ‘Q- 
operator’ method, is used to obtain an expression for the complete generating function of the cumulants 
of the current of the open ASEP, with the exact same structure as for the periodic case [32], proving a 
conjecture previously emitted in [127] . We will give a brief account of that method in section 


Many variants of the ASEP have also been studied. A matrix Ansatz, akin to the one we mentioned 
before, was found for the steady state of the periodic multispecies ASEP |128[ I129j . The case of a 
single defect particle was analysed, for itself |130H132j or used as a way to mark the position of a 
shock |133L 1134] . Different updates procedures were considered and compared for the discrete time case 
[13511136j . A system with two interacting chains was studied in |137j . The ASEP was also considered 
with entry and exit of particles in the bulk of the system [138] . disordered [139] or smoothly varying 
[140] jumping rates, a single slow bond [14111142) . repulsive nearest-neighbour interactions [143j . or on a 
two-dimensional grid [14411145] . 

Finally, on the numerical front, the ASEP has been used to develop and test numerical algorithms 
aimed at producing and analysing rare events, such as a variant of the ‘density matrix renormalisation 
group’ (DMRG) algorithm, first in [146j and later in [14711148] . as well as numerical implementations of 
the MET [14911151] . and the so-called ‘cloning algorithm’ [15211155] . 


IVB2 


B. Steady state of the open ASEP 

Before we look into the fluctuations of the current in the ASEP, let us first see what can be said of the 
average current and of the steady state of the model. We will start with a simple mean field approach, 
which gives good results in the large size limit, and we will then present the exact expressions of these 
quantities, in therms of the so-called ‘matrix Ansatz’ [33] . 


1. Mean field 


As we recall, the master equation reads: 


||P*)=M|P,) 

where M is the Markov matrix of the open ASEP: 


(52) 


L-1 

M = Too + Mi + rriL (53) 

i=l 


with 





0 

0 

0 01 



Too = 

—a 7 
a —7 

, M, = 

0 -g 1 0 

0 g -1 0 

, rriL = 

^7 

1 

1_ 




.0 0 

0 

0 




(54) 


We shall write the configurations of the system as C = where rii S {0,1} is the occupancy of 

site i. If we trace equation (52) over all n^-’s except for one at site i which is taken to be 1 (which means 
projecting it onto (l|(5„;q), we get an equation for the time evolution of the mean density at that site: 


dt 


(n,) = (l|<5„.,iM|Pt). 


(55) 


The term only affects two matrices from the sum in eq.(53l, and each of the matrices is stochastic 
individually, so that only a few terms remain in the right-hand side. A straightforward calculation yields 
the following expression: 


(^z) — Ji—1 Ji 


( 56 ) 
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with 


Jo = a((l - ni)) - 7 (ni), 

Ji = - rii)) - g((l - n,_i)ni), 

Jl = PiriL) - i5((l - ul)). 


(57) 

(58) 

(59) 


Equation (56 1 shows, as expected, that matter is conserved in the system: the variation of density at a 


site is equal to the current coming from the left minus the current leaving to the right. No approximation 
has been made so far, but we cannot solve these equations as they are: each current involves two-point 
correlations, so that the evolution of the densities {rii) is not autonomous. In principle, one would 
then have to compute, in a similar way, quantities such as -^{nirii+i), which would involve three-point 
correlations, and so on until one reaches the final L-point functions [10511107] . 

We will not be so ambitious, and will instead do a brutal approximation which will allow us to 
solve everything directly: we will assume that all the local densities are uncorrelated, which is to say 
(riinj) ~ {ni){nj). The resulting system of equations is autonomous. We are interested in the steady 
state, for which all time derivatives are zero, so that eq.(56) imposes that all the currents Ji are equal 


(in a one-dimensional system, a gradient-free current needs to be constant). Equations (57 
become 

J = a{l- (ni)) - 7 (ni) 

= (ni_i)(l - {rii)) - q{ni){l - {ui-i)) 

= P{nL) - i5(l - (rii)). 


591 then 


(60) 

(61) 

(62) 


The second of these equations gives a recursion relation between {rii) and (ni_i), which can then be used 
L — 1 times to express (n^) as a function of (n-i). Then, the first and last equations can be used to get a 
single equation on J, which fixes its value as a function of all the parameters of the system (L, q and the 
four boundary rates). These calculations can be found in |105j for the TASEP. We are only interested 
in the large size limit, so we will instead use a method similar to that of |106j . 


Writing {rii) = p C and taking L to infinity, equation (61) becomes 


J = (1 - 9)p(l - P)- 


with 


_p(^) = -V J = 0. 


(63) 


(64) 


Note that we have rescaled time in this last equation, by a factor L, and that we are keeping the gradient 
in eq. (63) even though its pre-factor goes to zero, as this will give us more information on the shape of 


the density profile for large sizes, and is in fact necessary to find the correct steady state. 


There remains the problem of finding the correct boundary conditions from eqs. (60) and (62), assuming 


that only the boundary rates from each side influence the corresponding boundary condition. This can 
be done by considering the situation where the steady state is homogeneous: p{x) = p. Equations (60) 
and (61) then give 


which is solved hy p = Pa = wrv, with 


1+a 
1 ■ 
' ^. 


a{l - p) - 7 p = (1 - q)p{l - p) 

1 

(1 — g — a -I- 7 ) -I- \J{1 — q — a + 7 )^ -I- 4 a 7 


Doing the same at the left boundary, we get a density pb = with: 

1 


b = 


2/3 L 


(1 - g - /3 -I- (5) -I- v^(wJgJw3Jr(3)JT4(^ 


(65) 


( 66 ) 


(67) 


Those two densities Pa and pb can be considered as the effective densities of the reservoirs to which the 
system is connected. We will take those densities as boundary conditions p are p(0) = pa and p(l) = pb 
in all future calculations. 
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Now that the equation and the boundary conditions are set, it is time to solve it. In principle, the 
standard way to do this would be to propagate the left boundary condition pa to the left side of the 
system through eq.(63l, keeping J as an unknown, and identifying the result with ph, thus obtaining an 


equation for J in terms of the boundary densities. Plugging this back into eq.(63l and solving it then 


yields p{x). It is in fact much simpler to fix J instead, to find all the density profiles compatible with 
that value, and to determine which boundary conditions are appropriate only at the end. 

Looking at equation (63), we see that the sign of Vp depends on the difference between J and (1 — 
q)p{l — p). Moreover, the gradient of p needs to be of order L to compensate for any finite difference 
between these terms. We can therefore argue that J cannot be larger than (which is the maximal 
value taken by (1 — 9)p(I — p)) or smaller than 0, otherwise |Vp| would be larger than some constant of 
order L, and p would diverge. This means that there is a density pc < 1/2 such that J = {l — q)pc{\- — pc)- 
We then have (fig.-[^: 


Vp < 0 

for 

P < Pc, 

(68) 

Vp > 0 

for 

Pc < p < (1 - Pc), 

(69) 

Vp < 0 

for 

(1 - Pc) < p. 

(70) 


which is to say that p gets away from pc and closer to (1 — pc). 



FIG. 5. Variations of p depending on its position with respect to pc and 1 — pc. For x increasing, (1 — pc) is an 
attractive fixed point, and pc is repulsive. 


Moreover, it is straightforward to check that p{x) is a hyperbolic tangent between pc and I — pc, and 
approximately an exponential otherwise, with a scale . On hg.j^ we draw some of the possible profiles 
for a certain value of J < . 

Since Pc is repulsive and 1 — Pc attractive, there are in fact only a few possibilities: 

• Pa = Pc and pb < 1 — Pc so that J = {1 — q)pa{l — Pa) ; this is represented in orange on fig.-[^ 
and requires Pa < \ and Pa < 1 — Pb ; this is called the Low Density phase: the left boundary 
imposes its (low) density to the whole system, apart from an exponential boundary layer at the 
right boundary. 

• Pb = 1 — Pc and Pa > Pc so that J = (1 — 9 )pb(l — Pb) ; this is represented in red on fig.-[^ 
and requires pb > | and pa > 1 — Pb ; this is called the High Density phase: the right boundary 
imposes its (high) density to the whole system, apart from an exponential boundary layer at the 
left boundary ; all this can be obtained from the Low Density phase through a leftoright and 
particleohole symmetry. 

• Pa = Pc and ^ = I — Pc so that J = (1 — q)pai^ — Pa) = (1 ~ Q)Pbi^ — Pb) ; this is represented in 
green on fig.j 6 and requires po = I — Pb < 5 ; this is called the Shock Line, since the domain wall 
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FIG. 6. Possible density profiles for a given pc- All the profiles with pa in the red region converge to p;, = 1 — pc. 
All the profiles with p(, in the orange region come from pa = Pc- 


going from pc to 1 — pc is a shock, which can be positioned anywhere in the system ; the steady 
state is in fact a superposition of all possible shock positions, and the average density is linear from 
Pa to Pfo. 

• The only situation which is not accounted for by these three cases is pa > | and pb < \ ] since p 
cannot decrease between pc and 1 — Pc, this requires Pc = 1 — Pc = 51 so that J = ; this being 

the largest current allowed, this is called the Maximal Current phase ; the density is close to ^ in 
the whole system, apart from algebraic boundary layers on both sides. 

We can summarise this by drawing the phase diagram of the system (fig.-[^. The transitions between 
the MC phase and the HD and LD phases are continuous in both the current and the density profiles. 
The transition over the SL, however, is discontinuous in the profiles (the mean density goes from pc to 
1 — Pc), but still continuous in the current. 


2. Matrix Ansatz 

In this section, we present the exact form of the steady state which we approximated in the previous 
one. It is formulated in terms of the famous matrix Ansatz, devised by Derrida, Evans, Hakim and 
Pasquier in |39j . using the recursion relations found by Derrida, Domany and Mukamel in |105j for the 
TASEP. Since the techniques involved here are rather specific to the ASEP, as they are related to it being 
integrable (although matrix Ansatze have been used for models which are not known to be integrable, 
such as the ABC model |I56j l. we will be much less precise here, and will merely give the main results 
as well as a few indications on the methods appropriate to approach them. For more details, one may 
refer to [31] or to section H. 2.1 of [T|. 

The statement of the matrix Ansatz is the following: for any configuration C = {rii}, with rii € {0,1}, 
the steady state probability P*{C) can be written in terms of a product of L matrices which may take two 
values D and E, sandwiched between two vectors ((IT|| and ||P)), up to a normalisation (those vectors 
are written using double bras and kets to distinguish the inner space of D and E from the physical space 
on which M acts). The product of matrices corresponding to C is obtained by multiplying matrices D 
for each particle, and E for each hole, in the same order as they appear in the configuration. In other 
terms: 


P^iC) = ^((W|| n {n,D + (I - n,)E) ||1/)) 


with the normalisation factor equal to the sum of all possible products, which is to say 


ZL = iW\\{D + E)^\V)) 


(71) 


(72) 
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FIG. 7. Phase diagram of the open ASEP. The values of the mean current in each phase are given in blue, and 
the mean density profiles are represented in the insets. 


SO that '^P*{C) 
c 

size 6 is given by 


= 1. For instance, the stationary probability of configuration 101101 for a system of 

{{W\\DEDDED\\V)) 

{{w\\(D+Er\\v)) ■ 


These matrices and vectors are of course not arbitrary, but must verify the following conditions: 


{{W\\{aE-jD) = {l-q){{Wl 
DE - q ED = {1 - q) {D + E), 
iPD-SE)\\V)) = il-q)\\V)). 


(73) 

(74) 

(75) 


Given these, it is rather straightforward to check that M\P*) = 0. 

We can obtain the stationary current from any of the equations (57 
corresponding equation from (73 - 73). For instance, equation ([5^ writes 


59), in combination with the 


J=^JW\\{D+Ey-^DE{D+E)^-^\\V))-q^{{W\\{D+Ey-^ED{D+Ef-^V)) = ( 1 -?)%^ ( 76 ) 


using equation (74) to simplify the expression. In all cases, we obtain 


J = {l-q) 


Zl-1 

Zl 


(77) 


and we are left having only to calculate for any L. 


This calculation was first done in [^, while the simpler equivalent for the TASEP can be 
[39]. Since Zl is a projection of {D + E)^, the natural procedure is to diagonalise {D + E). 
two new matrices d and e such that D — 1 + d and E = 1 + e, the bulk algebra (74) becomes 


found in 
Defining 


de — q ed = {1 — q) 


(78) 


which is that of a g-deformed harmonic oscillator where d is the annihilation operator, and e the 

creation operator. The eigenstates of {D + E) can then be found to be g-deformed coherent states, with 
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a unitary complex parameter z. One can then write the corresponding representation of the identity into 
Zl) and obtain, after a few lines of calculations, 


2 Jg 2iTTz (az, ajz^ dz, a/z, bz, bjz, bz, b/z)^ 


(79) 


where (•)oo is the g-Pochhammer symbol 


(a;)oo = 


(80) 


k=0 


with the notation convention [x,y)oo = {x)oo{y)oo- The two parameters d and b are defined similarly to 
a and b, but with a minus sign in front of the square roots, and have absolutely no importance in any of 
the calculations that we will do. 

For a < 1 and 6 < 1, the domain of integration is the unit circle, which is to say that we take the 
residues at all the poles of the integrand which are in S' = {0; aq^, dq^, bq^,bq^}ken, and not at the other 
ones (which are their inverses). Since is analytic in all the parameters, the poles of F at which we 
have to take a residue are always the same, even if one of them leaves the unit circle. For that reason, 
the integral in (79) must be done around S rather than the unit circle. 

We may rapidly remark on the origin of each term in that expression: (1 + z)^{l + z~^)^ is the 
eigenvalue of {D + E), to the power L ; (z^, z“^)oo comes from the normalisation of the eigenvectors of 
{D + E) ; (az,a/z,dz,d/z)ao and [bz,b/z,bz,b/z)oo come, respectively, from the scalar product between 
((fF|| and the right eigenvector of [D + E), and between ||P)) and the left eigenvector of {D + E). 


We can finally take L to infinity in Zl, to obtain the average current for large sizes. First of all, since 
the contour integral is done around an infinite number of poles, we will keep things as simple as possible: 
for any finite value of a and &, S contains all the poles inside the unit circle, plus a few poles outside of 
it (those for which aq^ > 1 or bq^ > 1), minus the inverse of those poles. Because of the symmetry of 
the integrand, the poles to not be taken in the unit circle have the opposite residue of those to be taken 
out of the unit circle, so that, all in all, the integral can be written around the unit circle, plus twice 
the residues around every pole of the form aq^ > 1 or bq^ > 1 outside of the circle. The poles related to 
d and b always stay inside of the unit circle, which is why they do not matter. This is summarised on 

fig-! 


l/a I 1/go 

® i o- 

q^a qa i 


1/g^a 


—o- 




FIG. 8. Only the pole at 0 and those related to a are represented ; the dashed line in the hrst figure is part of 
the unit circle ; the three sets of contours are equivalent. 


To find the large L behaviour of the integral, we need to note that, be it on the unit circle, or on the 
line [1, +oo[ on which all the poles sit, the dominating part of the integrand is (1 + z)^{l + z~^)^ which 
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is largest on the real axis and increases with z. If any of the poles in S are outside of the unit circle, the 
one which is the furthest from 1 dominates, so that Zl ^ max[(l + a)'^(l + (1 + 6)^(1 + 

Otherwise, the integral on the unit circle is dominated by the value at z = 1 and Z]^ ~ 4-^. 

It is left to the reader, as an exercise, to replace these expressions in eq. (0 and recover the phase 
diagram of J. 

Using similar techniques, one can compute the average density, as well as correlation functions |64l 
fTIMfTTTII . and validate the mean field calculations performed earlier. Note that we will treat contour 
integrals of the kind we just encountered in more detail in the next section, where we will see that all 
the cumulants of the current involve integrals of powers of the very same integrand. 
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IV. EXACT CUMULANTS AND LARGE SIZE LIMIT 


In the previous section, we saw what could be said of the average current in the open ASEP. In this 
section, we use the methods presented in m to extend our analysis to the fluctuations of the current 
as well, through the computation of the generating function of its cumulants, which, as we saw, is the 
largest eigenvalue of a deformed Markov matrix. We will first construct that matrix, and remark on 
some of its symmetries, which will help simplify the problem. We will then see how to obtain exact 
expressions for its largest eigenvalue, using techniques relying on the integrability of the model: this will 
be done for the periodic TASEP, using the coordinate Bethe Ansatz [S5], and for the open ASEP, using 
the Q-operator method at only a moderate level of detail, as these methods are entirely specific to 
the ASEP and are useful for other integrable models but not for other driven lattice gases. Once these 
exact expressions have been obtained, we will see how we can take the limit of large sizes and obtain the 
behaviour of the large deviation function of the current for small fluctuations. 


A. Deformed Markov matrix for the current 


The very first thing we need to define is precisely which current we intend to analyse. As we saw in 
section |IIIB[ we can define a current for each bond in the system. Their averages Ji need to be equal 
in the steady state, but we can in principle choose which current to monitor, take any linear combination 
of them, or even keep them all as independent variables. 

Each current defines a time-additive jump observable with U{C',C) = ±1, in the notations of eq.(20), 
if the transition increases or decreases that specific current, or 0 if the jump is at a different location. We 
can then deform our Markov matrix with respect to all these currents, each with a conjugate parameter 
/r, (fig.-[§. 


Jo 


a 


h 


-Mo 


■7 




JL = /?. 




-ML 


A 


n 




-fJ-i 


ji ~ 1^? I Q I 


FIG. 9. 


We obtain the following matrix 


L-l 

= »^o(mo) + X! + "^l(m/) 
2=1 


with 


wo(mo) = 




ro 

0 

0 

O' 



—a je 
ae^° —7 

, Mi{pLi) = 

0 -q 

0 ge-'^- 

e^. 

-1 

0 

0 

, mLifJ-L) = 

—6 j5e^^ 

6e~^^ —/3 



0 

0 

0 

0. 




(81) 


(82) 


(where, as before, it is implied that mo acts as written on site 0 in the basis {0,1} and as the identity 
on the other sites, and the same goes for on site L\ similarly. Mi is expressed by its action on sites i 
and i-|-l in the basis {00, 01,10,11} and acts as the identity on the rest of the system). 

The largest eigenvalue E{\pL^) of that matrix is the joint generating function of the cumulants of all 
the local currents, and the left and right eigenvectors carry the probabilities of configurations conditioned 
on the values of the integrated currents going to or coming from the steady state (as explained in section 
0. 
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Before proceeding to analyse that deformed Markov matrix, we should remark on a rather useful 
symmetry. If one considers the diagonal matrix i?i(A) (with ^ < i < L) with an entry e^ for all 
configurations for which site i is occupied, and 1 otherwise, one may easily check that the matrix 
similarity simply replaces and Mi{iJ,i) by, respectively, — A) 

and + A), and leaves the rest of unchanged. That is to say that part of the deformation 

is transferred from to Using combinations of this transformations for any sites and 

parameters A, we conclude that all the Markov matrices deformed with respect to the currents are similar, 
and therefore have the same eigenvalues^ as long as the sum of the deformation parameters p. = J2i=o h'i 
is fixed. Note that the eigenvectors, however, are different, but related to each other through those simple 
transformations. 

We will therefore write E{p) instead of E{{pL^). Since we are mainly interested in the eigenvalues 
of rather than its eigenvectors, we may choose a specihc combination of /i^’s to simplify our 

calculations. We will therefore put all the deformation on the first jump matrix mo(/i) and leave the 
others un-deformed, unless specified otherwise. We will write 

L-l 

= mo(Ai) + ^ + ruL- (83) 

2=1 


We may make one last remark on the deformed Markov matrix. There is a particular set of weights 
{/ii} dehned by 


a 


{fj,o = i^logl - fJ.i = V log {I/q), fj.L = ly log [-]} 


for which becomes: 


Too = 




o 

o 

0 01 



-a 

-7 

, Ml = 

0 -q q-^ 0 
0 -I 0 

, rriL = 

-5 

_p 



O 

o 

o 

o 




(84) 


(85) 


which is the deformed Markov matrix measuring the entropy production. We see immediately, as before, 
that 




( 86 ) 


which implies the Gallavotti-Cohen symmetry for the eigenvalues and between the left and right eigen¬ 
vectors of Ml, with respect to the transformation n <->■ (—1 —i^). 

Considering that p, = vlog ) i we also obtain the Gallavotti-Cohen symmetry related to the 

current, namely 


U(/r)=u(^-log( 


a/3 


l^q 


L-l 




(87) 


which is also valid for the other eigenvalues of M^, and the corresponding relations between the right and 
left eigenvectors, as well as a simple relation between the microscopic entropy production s, conjugate 
to V, and the macroscopic current j, conjugate to p: 

ajd 


s= j 


log(, 


There are several points to be noted here. First of all, those weights are ill-defined for the TASEP: 
micro-reversibility (i.e. the fact that for any allowed transition, the reverse transition is also allowed) is 
essential to have a fluctuation theorem. Moreover, if we take either the q —0 or the L ^ oo limit, the 
centre of the Gallavotti-Cohen symmetry p = —A log ) is rejected to —oo, so that the ‘negative 

current’ part of the fluctuations is lost. 

Finally, we may consider the detailed balance case, where = 1- In that case, we saw in section 

IIB 4 that there is no entropy production whatsoever, i.e. that E^v) = 0 (where we use the letter E 


abusively, since it is not the same function as E{p)). We see, indeed, from eq.(88), that s = 0. This does 


not mean, however, that j = 0: the deformations through p and u are in that case not equivalent, and 
E{p) 0. The only implication this has on E{p) is that it is an even function: E{p) = E{—p), all the 
odd cumulants are zero, and positive and negative currents of the same amplitude are equiprobable. 
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B. Exact cumulants of the current from integrability 


As we have mentioned earlier, the deformed Markov matrix of which we seek to extract the largest 
eigenvalue is integrable: its algebraic properties make it tractable in a systematic way, at least in principle, 
much like a quadratic one-dimensional spin chain Hamiltonian can always be diagonalised through free 
fermion techniques, although diagonalising it in practice is still a difficult problem. A wide variety of 
methods have been developed to tackle integrable models, such as the Bethe Ansatz and its variants 
(coordinate |158) . algebraic [55], functional [35], off-diagonal [15911160] . modified [161j b the Q-operator 
approach [53] , the q-Onsager algebra approach [16211163| , the separation of variables approach [164| , and 
more [US]. 

In this section, we will focus on two of these approaches, which are well adapted to our endeavour. We 
will first examine the coordinate Bethe Ansatz for the periodic TASEP [SS] , which is one of the simplest 
methods among the ones we mentioned, but requires the number of particles in the system to be fixed. 
We will then look at the Q-operator approach for the open ASEP [126| . which allows to calculate the 
generating function of the cumulants of the current exactly, although the analytical form of the solution 
has a few drawbacks, as we shall see. 

Considering that these methods and calculations are rather specific to the ASEP and its being inte¬ 
grable, and are of little use for other bulk-driven particle gases, we will only give the minimum level of 
detail necessary to understand how the results are attained and why they have their peculiar structure. 
References will be given along the way for the readers in need of more detail. 


1. Periodic case: coordinate Bethe Ansatz 

The coordinate Bethe Ansatz is perhaps the simplest way to approach the ASEP [15811166] . but requires 
the number of particles to be fixed. It relies on the fact that integrable systems can be understood 
as a generalisation of free fermions, where the anti-commutation rule depends on the particles being 
exchanged. As a result, the eigenstates of the system can be written as generalised determinants of 
one-particle eigenstates, where the coefficient of each term is a function of the permutation instead of its 
sign. 

We will see here how this can be used for the totally asymmetric case, as was first done in [96], and 
which, while being much simpler than the general case, has essentially the same behaviour. These results 
were later extended to the general periodic case in [55] . 

To make our calculations simpler, we will count the current on every bond in the system with equal 
weight, so that the system retains its translational invariance. The deformed Markov matrix we are 
considering is therefore 


with 


L 

2=1 




0 0 0 0 

0 0 0 

00 - 10 ’ 

0 0 0 0 


(89) 


(90) 


acting on sites i and i + 1, with L -|- 1 = 1. Note that the number of particles, N, is conserved. 

Let us first take A^ = 1. The system is then simply a totally biased random walk on a circle. Its 
eigenvectors are of the form 

|^«(z)) = ^z"|x), (91) 

X — 1 

where \x) is the configuration with the particle at site x. The periodicity condition imposes that = 1, 
which has L solutions. The eigenvalue associated with each of these solutions is A(z) = -l). 
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We now consider N = 2. Around configurations where the two particles are far from each other, at 
positions xi < X 2 (with an arbitrary site being labelled as 1), the system looks like two independent 
random walks, so that states such as 


^ (92) 

Xi <X2 


are locally stable, with an eigenvalue A(zi, Z 2 ) = —h ^ -2). This cannot be an eigenstate, since 

the two particles are in fact not independent, so that configurations where the two particles are next 
to each other would not be multiplied by K{zi,Z 2 )'- the contributions that would have come from the 
configuration where both particles are at the same site are missing. However, we see that the eigenvalue 
is invariant under exchange of zi and Z 2 , so that there are (at least) two locally stable states with that 
same eigenvalue. Taking a suitable combination of the two allows to make those missing terms cancel, 
and we obtain a true eigenstate. This is somewhat similar to the ‘method of images’ which is generally 
used for random walks with walls. Let us therefore consider 

\'ip’^'^\zi,Z 2 )) = ^ a z'l^z^'"+b Z 2 ^zl^\xi,X 2 ). (93) 

X\<X2 


Writing the eigenvalue equation for a configuration + 1), the coefficients a and b need to be such 

that 

A(zi, Z2)(a zl^z^^+^+b z^^zl^+^) = e'‘/^(a z^^-^z^^+^+b (a z^^z^^+^+b z^^z^^+^) (94) 

which simplifies to 

- Z2) + b{e^/^ - zi)=Q. (95) 

Moreover, due to the periodicity of the system, and because the eigenstates of cannot depend on 
the arbitrary labelling of the sites which we chose, we need to have 


a + b zf = a z^^z^^+^ + b zf zf 


(96) 


for any xi and X 2 , which is to say 


a = b Zi and b = a z^. 


(97) 


Since we are only interested in the eigenvalues of M^, we do not need to worry about a and b. Eliminating 
them from eq.(97), we get the ‘Bethe equations’ 


= to, . = 


(e'"/-^ - Zj) 

of which the solution can then be used to obtain A(zi, Z 2 ). 


(98) 


This simple case can then be extended to more particles without much effort: the integrability of the 
model ensures that triple collisions (when three particles are on adjacent sites) are merely combinations 
of double collisions, and do not add extra constraints to the z^’s. It is then straightforward to generalise 
eq.(95) to a relation between coefficients of terms with only two adjacent Zi’s being exchanged. Elim¬ 


inating them from the relations imposed by periodicity, we obtain the Bethe equations for N particles 
(and we recall the expression of the eigenvalues): 



For more details on how to obtain these equations, one may refer to chapter 8 of [521 or to [T]. 


(99) 
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The Bethe equations are a system of N coupled non-linear equations, which cannot be systematically 
solved in principle. Moreover, since the number of solutions to these equations is unknown a priori, it is 
not obvious that every eigenstate of the model corresponds to one of these solutions. 

That being said, the fact that we are looking for the largest eigenvalue of M^, which has non-negative 
entries save for its diagonal, ensures that we can attain our goal: the Perron-Frobenius theorem 
tells us that the largest eigenvalue E{pL) of is always non-degenerate, which means that we can follow 
it, and the corresponding eigenvectors, continuously by varying /i. Moreover, we know that £’(0) = 0, 
and it is easy to check that the right eigenvector for that eigenvalue at /r = 0 is uniform in the chosen 
occupancy sector (i.e. all the weights of configurations with N particles are equal, and all the others are 
0) ; we will write it as IIat). It corresponds to a Bethe state with all z^’s equal to 1. This information, 
combined with the Bethe equations, is enough to obtain an expression for 

To that purpose, we first need to change variables, for a simple matter of convenience: we will write Zi = 

eA‘/i(l-P 

so that the eigenvalues and the Bethe equations become, after a few simple manipulations, 


N 


A({^J) = E 


1 + 2/i 


- 1 


with 


.a + y^y 


Vi 


N 


= (-i) 


N-1 


N 

n% 

■i=i 


-1 


for 


= l..iV. 


( 100 ) 


The state we are interested in then corresponds to the solution with —>■ 0 when p ^ 0. Notice that 

N 

the right-hand side of the Bethe equations is the same for every i. Let us define B = yj, 

i=i 

and 


h{y) = 


(1 + 2 /)' 

o,N 


( 101 ) 


The Bethe equations then tell us that all the yi’s are roots of 1 — Bh{y), with the self-consistency 
condition given by the definition of B. Moreover, it is easy to check that, for y small enough, and B 
smaller than 1,1 — Bh{y) has exactly N roots inside of the unit circle among its L roots. Since we know 
that the y^’s should be small, these are the roots we are looking for. 

Summarising all this, we have three relations which we need to combine in order to obtain an expression 
of E{fj,): 

• the Pi’s are the roots of 1 — Bh{y) which are inside of the unit circle ; 

N , , 

•^^(m) = E(t4-i) ; 


N 


fi= log(-£«/y*). 


N 

Since both E and ji are of the form ^ we can write them as contour integrals around the unit 

i—1 

circle ci, using the first relation: 


N 




2 = 1 


dz {-Bh'{z)) 

z27r 1 — Bh{z) 


( 102 ) 


which, after an integration by parts, becomes 

Ef^,)iog{l-Bh{z)) 


(103) 


with f'{z) = —(l-|-z)“^ for E and f'{z) = —z~^ for /r. We can also check that there is no extra constant 
part coming from the integration by parts of a logarithm: for /r —>■ 0, both E and y have to vanish, which 
is the case here through B —>■ 0. 

We therefore have an implicit expression for E{y) in terms of B, given by 




( 104 ) 
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and 




( 105 ) 


which we may expand as series in B, and calculate the coefficients (which are all binomial numbers ; it 
is left as an exercise to the reader to check it): 


and 


M = 


k=l 


Bk 

~k 


with Ck 



E{y:) 


k=l 


Bk 


with 


/kL-2\ 

\kN-l)' 


(106) 


(107) 


The periodic ASEP, with q ^ 0, yields results with a similar structure, as was found in 


From these formulae, one can, in principle, calculate any cumulant in a finite number of steps, by 
inverting qi{B) order by order and injecting it into E{B). Since we want to take a Legendre transform of 
i?(/i) to obtain the large deviation function of the current, we are rather interested in an closed expression 
of E{^), at least in some limit. We will see how to obtain that shortly, but we will first have a look at 
the open ASEP, where the coordinate Bethe Ansatz only applies in special cases [331 [Ml I3H] • 


2. Open case: Q-operator method 


We now go back to the open ASEP, with a deformation on mo as in eq.(83). Because the number of 
particles is not fixed, the coordinate Bethe Ansatz cannot be applied for generic boundary parameters, 
and we will have to use a different method: that of the so-called ‘Q-operator’, sometimes called ‘aux¬ 
iliary operator’ nsH], which gives a relatively simple way to obtain the expressions we seek. The main 
difference between this method and that of the Bethe Ansatz is that it does not require an Ansatz for the 
eigenvectors to access the eigenvalues. Moreover, it is entirely algebraic: it yields relations verified by 
the matrix rather than by parameters entering an Ansatz for its eigenvectors (which is the case for 
the Bethe equations), so that the completeness of the solution does not need to be proven a posteriori. 

The calculations and proofs involved in applying the Q-operator method are rather lengthy and tech¬ 
nical, and we will not dwell on them here. All that we will see here is based on UliMj, to which the 
reader may refer for more details. 


At the core of this method is a generalisation of the d and e matrices which we used in section |IIIB 2| 
to build the Matrix Ansatz. Let us consider similar matrices, with a slightly different q-commutation 
rule, as well as a third matrix A, defined as 

de — q ed = {1 — q){l — x^A^) , dA — qAd = 0 , Ae — qeA = 0, (108) 


where a; is a free parameter. 

As before, one can think of d and e as the annihilation and creation operators of a q-deformed harmonic 
oscillator with an extra parameter, and A as the q-counting matrix g" (where n is the number of 
excitations). In this usual representation, they can be written down as 


OO OO 

d= ^(l-g")||n-l))((n|| , e = ^(1 - a;2g")||n-t l))((n 

n—1 n=0 


A = Y,q-\\n)){{n 

n—0 


(109) 


In section [Til B 2[ we used these objects to build a vector, with 1 + d 
Here, we need them to construct a transfer matrix, in a similar way. 
with those objects as entries: 


for a particle, and 1 -I- e for a hole. 
Let us define a 2 x 2 matrix X(x) 


X{x) 


1 -I- xA e 
d 1-1- xA ’ 


( 110 ) 
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expressed in basis {0,1}. We also need to define a matrix such that 


dA^ — e ^ A^d = 0 


Af_,e - e ^ eAfj, = 0, 


( 111 ) 


which is to say = e with the same structure as A but with q replaced by e Moreover, just as 


for the matrix Ansatz in section III B 2 we need to define special vectors to be placed at the boundaries, 
containing the parameters from toq and In the present case, we need four vectors ||y)), ((W||, ||F)) 
and {{W\, such that 


[/l(d + I + xA) - (5(e + 1 + xA) - (1 - 9 )] ||y)) = 0, (112) 

{{W\ [a{e + 1 + xA) — 'y{d + 1 + xA) — (1 — q)] = 0, (US) 

[f3{d — 1 — xA) — d{e — 1 — xA) + (1 — q)xA\ ||y)) = 0, (114) 

((W|| [Q:(e — 1 — xA) — 7 (d — 1 — xA) + (1 — q)xA] = 0. (US) 


As one can easily check, the first two are generalisations of the ones which we used for the Matrix Ansatz 
(which we recover by taking x = 0). 

Using all those elements, we build two transfer matrices U^(x) and T^{x) (each corresponding to a 
pair of boundary vectors), with a structure very similar to that of the matrix Ansatz: instead of a 
product of matrices D and E for particles or holes, we use the elements of X{x) for transitions between 
configurations (i.e. Wj at site k if the initial configuration has occupancy j at site k, and the final one 
has occupancy j at site k). We also need to add a matrix at the place where we are counting the 
current, which is between the left border and the first site in this case (if we had more general /r^’s, we 
would need a matrix A^. for each of those). For instance, the entry of U^{x) from configurations 001010 
to the right to 010111 to the left is given by ((IU||A^(1 + xA) d e d {1 + xA) d||U)). 

We write these matrices in the following way 


L 


C/^(x) = ((IU||A^n^^'^(^)l^))> 

(116) 

L 

r^(x) = ((iu||A^n^^'^(^)l^))’ 

(117) 


i=l 


with an exponent (i) on each X which serves only to signify that it acts on site i. Note that those 
transfer matrices are in principle defined up to a non-trivial normalisation, which can be a function of 
X, and is implicitly included in the boundary vectors. That normalisation plays a part in particular in 
equation (121), which is not homogeneous in 17^ and T^, and we choose it so that eq.(122| holds. 


The main property of these matrices is that any product t/^(x)T)j(?/) commutes with any other product 
Ufj,{x')T^{y'), and with M^: 

[U^{x)T^{y),U^{x')T^iy')]=0 , [U^{x)T^{y), M^] = 0, (118) 


which can be shown using so-called ‘R-matrices’, as is customary for transfer matrix methods |28j . Note 
that the dependence in x of each of the matrices is not only in the diagonal entries of X and in the 
relations defining the boundary vectors, but also implicitly in d and e through eq. (1081. 

Those two transfer matrices do not commute with one-another, but we may use them to define two 
that do: P{x) and Q{x) defined by 


P(x) 


U^.{x) U^{0) 


1 -1 


Q(x) = (l-e-^)-i[/^(0)r^(x) 


(119) 


such that: 


U^{x)T^{y) = (1 - e-nP{x)Q{y). (120) 

The downside of these is that, unlike 17^ and T^, P is not defined constructively because it involves 
[C/^(0)]“^ (which exists for a generic y), and does not have the same product structure. 
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Apart from their commuting with (which means they have the same eigenvectors), the operators P 
and Q have another useful property: given certain constraints on their variables, their product decouples 
into a transfer matrix with a finite auxiliary space and a shifted version of the product. More precisely, 
for every positive integer A:, 


P{x)Q{l/q'^-^x) = t^'^^x) + e-'^'^^P{q^x)Q{q/x). 


( 121 ) 


The matrix (x) is the transfer matrix of an integrable vertex model with a /c-dimensional auxiliary 
space. For instance, t^^'^{x) is the transfer matrix of the six-vertex model [53] ■ Moreover, t^^'>{x) is a 
scalar matrix, which we can calculate, choosing the appropriate normalisation for and 


s F{x) = (i + ») (!+»'-*) 

(ax, a/x, ax, a/x, bx, b/x, bx, b/x)o 


( 122 ) 


where the identity operator on c onf iguration space is implicit. We already encountered that function 
inside of a contour integral in eq.(79). t^‘^'>(x) and F(x) are strongly connected to M^, through the usual 


relation between the six-vertex transfer matrix and the Hamiltonian of the spin-^ XXZ chain (which is 
equivalent to the ASEP): 




F(qx) 


Combining eq.( 1211 for fc = 1 and fc = 2, we can eliminate P and obtain the ‘T-Q equation’: 

t’''^''(x)Q{l/x) = F{x)Q(l/qx) + e~^^F{qx)Q(q/x) 
which allows us to express in terms of Q alone: 


M.= i(l-,)dlog 


/ Q(q/x)\ 


\Q{l/x)) 

X — —1 


(123) 


(124) 


(125) 


This relation between and Q (which applies to the whole spectrum of A/^), along with eq.(121) for 
A: = 1: 


P(x)Q(l/x) = F(x)+e '^^P(qx)Q(q/x) 


(126) 


is enough to obtain E(fj,): they are a more complex version of the relations between E, fi and the Bethe 
roots we used in the periodic case (in fact, applying the Q-operator method to the periodic case, we find 
that the Bethe roots are the inverses of the roots of the eigenvalues of Q). 


From this point on, we need to use the ‘functional Bethe Ansatz’, as was done in |32j for the periodic 
ASEP, in order to unravel these relations and obtain an expression of E(^). We first need to define a 
matrix B which plays the same role as the intermediate variable we used in the previous section: 

B = -e2/^(Q(0))-^ = -e"'^(l - e-^){U^(0)T^i0)y\ (127) 


Moreover, just as in the previous section, we need to know the behaviour of B and Q in the /i —>■ 0 
limit. We find, from analytical calculations which we will not go into here but can be found in [T], that 
the first eigenvalue of B goes to 0, while the others remain finite, and that the roots of the first eigenvalue 
of Q(\/x) are inside of the unit circle ci if o < 1 and 6 < 1 (as was needed in section III B 21, while 
those of P are outside of it. Restricting ourselves to that first eigenspace from now on, this allows us to 
separate P and Q in eq.( 126) using a contour integral, and express in eq.( 1251 only in terms of F{x). 


In the following, the notations F, P, Q, and B will refer to the eigenvalues of the operators in the 
dominant eigenspace, and not to the operators themselves. 

Let us define a function W(x) as: 


W{x) 



P{x)Q{l/x) \ 


e~'^f^P(qx)Q{q/x) J ’ 


(128) 
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and a convolution kernel K, as: 


OO ^ 

^) = 2 E 




along with the associated convolution operator X: 


xim = 


dz 

^27rz 


f{z)K{z,z). 


(129) 


(130) 


Expanding log(P(x)) and log((5(l/a:)) in powers of x and 1/x respectively, we can easily check that 

-\og{P{qx)Q(q/x)/Q{G)) = X[W]{x), (131) 

which allows us to rewrite eq. ( fl^ as a functional equation of only one unknown function W: 


W{x) = ln(l - . 


(132) 


The last step is to take eq.(125) in the first eigenspace of M^, and eq.(128) at a; = 0, to find: 

dz 


E{y.) = \{l-q)j^\0g 


Q{q/x) 

Q{l/x) 


= hi-g) 


log 


z27r(l + z)2 \Q(l/z) 


Qid/z) 


fj, = -W{0) (133) 


where the contour integral form can be checked formally by expanding log{Q{l/x)) in powers of 1/a;. 
Considering that, for fj, small enough, P is holomorphic inside of the unit circle, we can replace 

I by —W{z) when expressing E{fj.) as a contour integral (since P will not contribute), and 


2 

obtain: 



dz 

i2'kz 


W{z) 


(134) 


and 




(135) 


All this was done for a < 1 and 6 < 1, but can then be generalised to any a and b through the same 
reasoning as in section |III B 2| for the mean current, replacing the unit circle ci by small contours around 
S = {0, g'=a, q'^d, q'^b, q^. 


As we see, the form of E(fi) is similar to what we found for the periodic TASEP. It is even more similar 
to the periodic ASEP case [32]: the expressions only differ by the factors 2 in AT and ^ in W{x), and by 
the function P{x) which is exchanged with h{x) as defined in eq.(101). 

Moreover, if we choose a = 1, d = —q, b = ^/q and b = —y/q, which is to say Q;=t, 7 =|, /3=1 
and 6 = q, special cancellations occur in F{x) which reduces to (1 + . This is the 

same as the function h for the periodic ASEP with 2L + 2 sites and L + 1 particles. Because of the extra 
factors 2 and t, the generating function of the cumulants of the current is half that which we found in 


the periodic case, taken at 2/r. This also works if we exchange a with b and d with b. Those two special 
points correspond to pa = ^ and 1 — pb = or the opposite, and are on the transition lines between 
the MC phase and the LD or HD phase. We will come back to this remark later. 


The expressions we just obtained are exact, and valid for any values of the parameters of the system, 
including its size. They are, however, somewhat unwieldy, especially if we want to perform a Legendre 
transform in order to obtain the large deviation function of the current, which is our goal. In the next 
section, we take the limit of large sizes and see how we can obtain closed expressions for E(p) and g{j) 
for small fluctuations of the current. 
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C. Large size limit 


In this section, we will need to take the result we just obtained for the generating function of the 
cumulants of the current in the open ASEP, and extract its large size behaviour, through various ap¬ 
proximations. We will not go into every fine detail of the necessary calculations, but we will endeavour 
to make them as easy to follow as possible. That being said, the calculations themselves might not be of 
interest to every reader, and we will make the results clearly visible at the end of each sub-section (i.e. 
they will be boxed). 


Before taking any limit, we will need to examine eqs.( 134) and (135) a little closer. The function W(z) 
that they contain is defined, in eq.(132), through a self-consistency equation. Expanding the logarithm 
in powers of B, we can express every coefficient by calculating W(x) perturbatively in B. The coefficient 
of B'^ will be a combination of k functions F, either at the same point or convolved together through K. 
For instance, the coefficient of B^ in fjL is 


1 

4 


dz 

i2ttz 


F{zf 


1 

4 


dzi 

%2 'kzi 


dZ2 

i2'KZ2 


F{zi)F{z2)K{zi,Z2) 


(136) 


As was painstakingly verified in it turns out that, in the large size limit, the term containing no 
convolutions is dominant in every case, and that, when all is said and done, the behaviour of E{^) differs 
from the TASEP case (where q — 0 so that K — 0) only by a global factor (1 — q). We will therefore 
save ourselves some trouble here and only consider that simpler case: q = j = 6 = 0. 


The simplified expressions we will have to examine are the following: 


^ TDk 

k=l 


with 


Ck = - 

k 2 


{0,a,fe} 


dz 

i2ttz 


F\zl 


(137) 




Bk 


with Dh = 


k=l 


dz 


2 J{o,aM *27r(l -H z)2 




where 


F{x) = 


(1-I-a;)-^)!-|-a; ^)^(1 — a:^)(l — a; ^) 


(138) 


(139) 


(1 — aa:)(l — a/a;)(l — 6a:)(l — b/x) 

As for the calculation of the mean current that we saw in section |IIIB 2[ the asymptotic form of these 
contour integrals, for L large, will depend on the position of a and b with respect to the unit circle (see 

Moreover, the correct way to handle that large size limit is not entirely straightforward, due to i?(/i) 
being expressed implicitly: each cumulant is a combination of C^s and D^s of same or lower order, such 
as 


E^t — 


iDiC^ - 2D1C1C3. - iD2CiC2 + 2D3CI 


(140) 


but the leading order of E^ in L might not be obtained by taking the leading orders of each term, as this 
may simply give 0 due to certain cancellations. As it turns out, this is in fact always the case: at leading 
order in L, in every phase, one finds Dk ~ JCk, where J is the average current, so that E(/r) ~ J/i. This 
cancels out in every Ej. other than Ei = J. For that reason, we will sometimes be calculating equivalents 
of 


E{^x) - Jfi = 


Bk 


with Dk = Dk — JCk 


(141) 


k=l 


instead, which we will find to be sufficient in every case we will examine. We will only ignore one point 
of the phase diagram: a = /3 = which is to say a = 6 = 1, at which the k first orders of Ek vanish. 

Finally, one should note that, because we will make approximations on the coefficients of E{p) ex¬ 
panded (implicitly) around fi — 0, the expressions we will obtain for E{fi) or g[j) will in principle be 
valid only for a small fz and a small fluctuation (j — J). 
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FIG. 10. Positions of a (blue circles) and b (red circles) with respect to the unit circle (black circles) in each 
phase of the open ASEP. 


1. Low/high density phases 


We start with the low density phase: a > 1 and a > b. The corresponding results for the high density 
phase are obtained by excha nging a and b. 

As we recall from section IIIB 2[ we can replace the integration set {0, a,5} by the unit circle plus 
twice the poles which are outside of it. The dominant part of each contour integral is then the residue 
around a, which is the one where F{x) is maximal. That residue is of order k in both Ck and Dk- We 
can therefore write 


OO 

m = -E 

k^l 

OO 

Eiti) = -Y, 


kl dz^~^ 1 Z i\z=a 

(142) 

fife ( />^-(^) 11 

k\ dz^~^ ^ (1 + z)^ / lz=a’ 

(143) 


with 




(z 


a)F{z) 


(1 + z)^{l + z-i)^(l - z^){l - z-^) 

(1 — az)(l — bz)(l — bjz) 


(144) 


One may recognise there a structure related to the Lagrange inversion formula [169] . Considering two 
variables w and z such that 


w = z + B(p{w) 

we can express a function / taken at w by expanding it around z, as: 


^ TDk f^k—1 

/(«■)-= 


/c=l 




Choosing /(z) = — log(z) for fj, and 1/(1 + z) for E{pl), we can then write: 

M = 

E{pi) = 


M =-log(w)+ log(a), 

1 1 


w + 1 a + 1 ’ 


(145) 


(146) 
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where w is as defined in (145). Combining those two equations, we get a closed expression for E{ii): 


E{fi) 


a e^ - 1 
a + 1 e^ + a 


(149) 


It only remains for us to perform a Legendre transform to obtain the large deviation function of the 
current: 


gU) 


{l-q) Pa-r + r(l - r) log(^ 


1 - /Oa 


Pa 



(150) 


where r is such that j = (1 — q)r{l — r) and we recall that pa = (iqyi) ■ 

This expression was first obtained in [88] using a method which we will discuss in section VI The 
formula (149), with an extra factor (1 — g) for the general ASEP, was obtained in |124j through a 
numerical approach to the Bethe equations. Contrary to what we noted earlier, these expressions are in 


fact valid in their whole phase, and not only for small p, as we will see in section VI 


2. Shock line 

We now look at the shock line, where a = 6 > 1. Unfortunately, the method we used in the previous 
section does not work here: the residue at a in Ck and Dk is of order 2k, so that we are missing all the 
derivatives of even order in p and E. Instead, we simply calculate the leading order of each coefficient. 
After simplifying every term as much as possible (see chapter IV of [T] for a detailed calculation), we 
obtain the large-size equivalent of our two series: 


2 a -b 1 ^ 


E{p)-{l-q) 


(1 + a)' 


rM = -(1 - g) 


La-1^^ (2k)\ 

9 „ “ pfe-2 

a?-I (2k)\ 


(151) 

(152) 


and we find that the cumulants of the current behave as 

for k> 2. 

These two series can be expressed in terms of the Lambert W function |170j . defined as the solution 
to a; = Wce^^. The series expansion of Wc{x) around 0 is: 


Wcix) = -J2 


-i-xY 


n—1 


so that p and E become: 


E{p) - (1 - g) 


P = 


2 Q. “h 1 


(l + a)' 


-p= {l-q) 


La- 1 1 
2 a 


yVc{VBl2) + Wc{-'/b/2) 


L2 a2 - 1 [2 (w£(^/2) + Wci-y^/2) 
+Wc{Vb/2)'^ + Wc{-VbI2)^ 


(154) 


(155) 

(156) 


Many things are known about W^, including its asymptotic behaviour, which is what we need. The 
main branch of the function, called Wo, is defined on the whole complex plane except for ] —oo,—1/e], 
and behaves as log(a;) (even for x complex, in which case the angular part of x can be neglected and we 
get log(|a;|)). This will be appropriate for B < 0. For L > 0, however, the functions yVc{—VB/2) in 
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FIG. 11. Plot of the Lambert W function. The principal branch (blue) behaves as log(a;) for x —>■ oo. The second 
branch (red) behaves as log(—x) for x —^ 0“. 


our expressions have to be continued analytically to the second branch W-i, on which x goes back from 
— 1/e to 0 and behaves as log(—x) (hg. 0- 

For B —>■ —oo, we therefore have >V£(±-\/i3/2) ~ log(|i?|)/2, so that: 


,^|^log(|B|) 

- (1 - q) , ° ~ (1 - log(|B|)V2 

(1 + a)^ — 1 


(157) 

(158) 


meaning that, for /i positive and small (since we are working in the L ^ oo limit ; see comment below), 
E behaves as 


E+{^i) - (1 - g) 


a 

(l + a)2 


/i~ (1-g) 


aia — 1) 9 

4()/TTF^- 


(159) 


We take the Legendre transform of this and get, for j > J = {1 — q)pa{l — Pa)'- 


9+ij) 


{J-Jf 

J(1 - 2pa) 


(160) 


(where J = (1 — q)pai^ — Pa) is the average current). 


careful analysis. Starting from equations (1371 and (138), which are exact for any L and B, we first took 
a large L limit and kept only the leading contributions, yielding equations (151) and (152), still valid 


The a priori domains of validity of equations (159) and (1601 are not entirely obvious, and require a 


for any B. We then took a large (negative) B, and obtained equations (157) and (158). Becau se of the 
order in which we took the limits, L is taken to be large enough so that the first correction to (151) be 


negligible for any B we consider, which is to say that L goes to oo faster than log(|i?|). For that reason, 
the equations we obtain here, and similarly in all the following calculations, are in principle valid only 
for p small, which is to say j close to its average. We will indeed see in section [VT| that these expressions 
correspond only to the leading order of those we find from the macroscopic fluctuation theory. 


To do the same for /j, < 0, we need to be on the second branch of Wc, for which we have B —>■ O’*". 
In that case, we have 1‘E) \og[B)/2, but Wci'/B/2) ~ 0 (that part being still on the main 
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branch). This gives us: 


(l + a)2 

so that, for /i negative and small, we have 


L a — 1 

- (1 - - (1 - log(|B|)V4 


f \ / A \ ^ {A \ f) 2 

- (1 - g) .. , ~ (1 - g L/ ,1X3 ^ ■ 


(161) 

(162) 

(163) 


We take the Legendre transform of this and get, for j < J: 


9-U) 


{J-J? 

2J(l-2p,)- 


(164) 


Notice that the dependence in L has vanished from both cases, so that even though all the cumulants 
of the current depend on L at ^ = 0, none of them do for a finite /r. Notice also that differs from 
g+ by a factor f, which comes from the fact that for /r > 0, both functions Wc in 9 and >V| in E{g,) 
contribute to the limit, whereas for p, < 0, only one of each does. This results in g{j) not being analytic 
at j = J, which is the signature of a non-equilibrium phase transition (see section 
of the phase on each side of the transition). 


VI for a description 


3. LD/MC transition line 


We now consider the transition line between the low density and the maximal current phase, for which 
0=1 and 6 < 1. As we noted earlier, the cumulants of the current on this line are related to those for a 
half-filled periodic TASEP of size 2L + 2. The calculations we will need to perform here are therefore the 
same as can be found in which we will reproduce, and use as a reference for the maximal current 
phase in the next section (which is similar but not identical). 

Taking the limit of large L in g and E (see chapter IV of [T] for a detailed calculation), we obtain 

L-i/2 - 

9 = 


2a/tt 


k=l 


1 r-3/2 °° 

E[g) - (1 - q)-,9 = -(1 - g)TWTV XI 


IG-^TT 




so that the cumulants behave as 


Ek ~ 7r(7rL)('=-3)/2 


Bk 

p/2’ 

(165) 

Bk 

p/2’ 

(166) 


(167) 


for k >2. Notice that the pre-factor L in (166) is a sign of the dynamical scaling 2 = | of the KPZ 
universality class. 


These series can be written in terms of the poly logarithm Li 5 / 2 (B), defined as 




k=l 


Bk 


r+°° , 

/ de e‘^\og[l-Be-^ ], 

J —oo 


(168) 


so that 


L-1/2 

E{g) - 

4 ^ 16V7r 


B H'{B) 

(169) 

H{B) 

(170) 
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As in the previous section, the cases /x < 0 and /x > 0 require different approaches. 


For /X > 0, we need to take B —)■ —oo. In this case, the integrand in H{B) can be approximated by: 

\og[l-Be-^"]^\og[\B\ e-®'] < log(|B|)] (171) 


where I[A] is the indicator of X, equal to 1 if A is true, and 0 if A is false. 

We can then estimate: 

H{B) ~ ^ de 0^(log(|S|) - e^) = ^ log(|i?|)5/2 


and 


so that 


B H'{B) 


Si/tt 


log(|B|)3/2 




1 /3\2/3_ 


-2/3 2j3,,5/3 


(172) 


(173) 


(174) 


We can then take the Legendre transform of this result. We find, for j > J = 


4 ■ 


9+U) ~ (j - 


57r(l — q)3/2 


(175) 


and notice that, for once, it depends on L. 


For /X < 0, we have to take into account the fact that the polylogarithm has a branch cut at 

X = 1 and is not defined for a; > 1. To remedy this, we have to go back to the expressions of /x and E 
in terms of the roots of (1 — Bh{x)) we saw in section IVB 1[ These expressions apply in principle to 
all eigenvalues of M^, and depend on which roots are included in the integral. In the case of the steady 
state, for B small enough, all the roots we consider are inside of the unit circle. However, it can be 
shown that, as B gets closer to 1, one of the roots zq goes to 1, and merges with its counterpart Zq^ from 


outside of the unit circle (fig.-12). Since we know, from the Perron-Frobenius theorem, that E(fj,) never 


crosses any other eigenvalue of M^, the choice of roots that consists in taking Zq^ instead of Zq must 
correspond to E{fi) as well (because they coincide for i? = 1). We can therefore find the correct analytic 
continuations for /x and E{p) in terms of B by finding zo, replacing its contribution in those series by 


that of Zq and taking B back from 1 to 0. This procedure is explained in more detail in 


We can find those two roots using equation (168). An integration by parts turns it into: 


H{B) 


J- 




Be-“ - 1 


(176) 


For 0 < H < 1, the poles in this expression are at 9± = EisJ— log(i3). The corresponding residues are 
(as explained in |171j L and we must subtract the one corresponding to and add the other 
one to a. We get, for B M 0: 


H{B) = ^v^[-log(H)] 


3/2 




BH'{B) = -4^71 [- log(H)]^ 
Putting those together, we find, for /x < 0: 


k=l 


fc3/2 


-4V7r[- log(H)] 


1/2 


- (1 - 


(177) 

(178) 

(179) 
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FIG. 12. Bethe roots for a periodic system with 20 sites and 10 particles. The roots are at the centres of the 
white discs. The unit circle is represented in black. On the left, where B < 1, a pair of roots can be seen to 
approach the unit circle on the real axis. On the right, where B — 1, those roots have merged. 


and, for j < J = , 


9-{j) 




8 

3(1 - 9)i/2- 


(180) 


Once more, we find a result which is independent of L. Moreover, we see that the phase transition 
that takes place here at /r = 0 is of a different nature than the one on the shock line: the behaviour of 
g{j) with respect to L changes from one side to the other. 


4- Maximal current phase 


Finally, we inspect the maximal current phase, for which a < 1 and b < 1. Once more, we start by 
calculating the large L behaviour of every Ck and Dk'- 


Eig) - (1 - 


9 = - 
-(1-g) 


£-1/2 

2i/7r 

£-3/2 

16v^ 


OO 


E 


OO 


E 


(2fc)! fc 

k\k(k+3/2) 

(2fc)! fe 

k\k(k+5/2) > 


(181) 

(182) 


so that, as in the previous section, 

iffe ^ (1 - g)7r(^£)('=-3)/2 


(183) 


for k > 2. 

We will use the same method as for the LD/MC line. We may still write 


E{f,) 


A* = 


£-1/2 

2i/7r 


B H'{B) 


I — q £ 3/2 
4 ^ 16y/^ 


H{B) 


(184) 

(185) 
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but this time, H{B) is defined as 




fc=l 


For /X > 0, i.e. B —)■ —cso, we have: 

log[l-B6»2e-®']~log[|B|6»2e-®'] l[\B\e‘^e-^^ > l]. 


( 186 ) 


(187) 


The upper bound of the integral can therefore be set at 9b such that \B\9^e~^‘^ = 1, in which we 
recognise the square root of the Lambert W function: 0b = yjB). For large B, it behaves 
as log(|B|)i/2, as we saw in section 


IVC2 


After estimating H{B) in the same way as in the previous 
section we find the exact same expressions for /i and A, leading to 


EM - ~ (1 - 


1 / 3 \ 2 / 3 ^ 


(188) 


and, for j > J = 


4 ’ 


9+U) ~ (j - 


32y/3L 
57r(l — g)3/2 


(189) 


This tells us that the phase which is reached by selecting positive fluctuations of the current is the same 
whether one starts from the inside of the MC phase or from its boundary. 


For /r < 0, the situation is slightly different from that of the previous section. The roots of (1 — BF{x)) 
behave similarly, but this time, there are two pairs of roots crossing the unit circle instead of one (fig.-131, 
both close to the real axis. Using the same procedure as before, we find them to have exactly the same 
behaviour with respect to S, the only difference being that we have twice as many residues as we had 
then. 



FIG. 13. Bethe roots for an open system with 9 sites. The unit circle is represented in black. This time, there 
are two pairs of roots that merge for a critical value of B. Those roots get closer to the real axis as L becomes 
large. 


This gives us: 

H{B) ^ ^y^[-log{B)f^ 


(190) 
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and 


BH'{B) - -8v^[- log(S)]^^^ 


so that, for /r < 0, 


E-(p) - ~ -(1 - 


and, for j < J = , 


5-(j) 




16 

3(1-9)1/2- 


(191) 


(192) 


(193) 


Note that we do not find the same behaviour for negative fluctuations of the current inside of the MC 
phase and on its boundaries, which indicates a phase transition between those two regions for /r < 0. We 
will explore that transition, and all those we mentioned in the previous sections, in more detail in section 


VI But before that, we will see what may be said of E{p) and g{j) in the limit of large fluctuations. 
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V. ASYMPTOTIC LIMITS FOR EXTREME CURRENTS 


In the previous section, we obtained the behaviour of the generating function of the cumulants and of 
the large deviation function of the current for small fluctuations. In particular, we found evidence for a 
number of dynamical phase transitions, in the different forms these functions take for positive or negative 
fluctuations. In this section, we examine the opposite limit: that of extreme fluctuations. We will do this 
by taking /i to ±oo in the deformed Markov matrix M^, and diagonalising the limiting matrix directly, 
without use of integrability-related methods. 


There are a few caveats which need to be pointed out before we start calculating. First of all, as we saw 


in section IV A the current in the open ASEP is proportional to the entropy production, and satisfies 
the Gallavotti-Cohen symmetry, which means that the /i —>■ —oo limit, corresponding to a large negative 
current, can be deduces from the fi —>■ +oo limit, corresponding to a large positive current. However, if 
we consider the totally asymmetric case instead, the Gallavotti-Cohen symmetry is destroyed, and the 
/i —>■ —oo limit becomes that of a null current (simply because negative currents are strictly impossible). 
Whether the behaviour of the TASEP for /x —>■ — oo and that of the ASEP for j —)■ 0+ (which corresponds 
to ^ log(Q;/?/7(5q^“^)+)are the same is not entirely obvious, but we will see evidence for it in section 

VI In the meantime, we will focus on the TASEP when taking the ^ > —oo limit. As for the /x —> -l-oo 
limit, we will immediately see that considering the ASEP or the TASEP makes no difference whatsoever, 
and we will focus on the latter simply for the sake of consistency. 

As for dealing with these limits in a proper way, one needs to be careful. Nothing guarantees that the 
equivalent of will be diagonalisable, and it will in fact not be in general. Consider for instance the 
case where we monitor the current on the first bond, and take /i to infinity. We are left with the matrix 


~ mo(/x) 


0 0 
ae'^ 0 


(194) 


acting only on site 1, which is not diagonalisable. This comes from the fact that the eigenvectors do not 
have good limits in that basis: some entries, to which vanishing elements of apply, will diverge, which 
makes it inappropriate to neglect those elements. However, if we find a basis in which the equivalent of 
is diagonalisable, then the eigenelements of the limit are the limits of the eigenelements. Luckily, we 
already have a set of basis transformations at our disposal, which consist in changing the way in which 
we measure the current (see section IVA). It turns out that putting a weight fj^i = on each bond is 
the simplest distribution which will yield diagonalisable limits. 

In this section, we will therefore consider the deformed Markov matrix defined as 


L-l 


= moifJ'o) + ^ Mi{ijLi) + mLifJ-i), 


(195) 


where 


moifio) = 


—a 0 

ae^° 0 


, = 


0 0 
0 0 


0 


0 0-10 

0 0 0 0 


, mLini) = 


0 /3e^^^ 
0 -B 


(196) 


with for all i’s, as our starting point. 


A. Low current limit 

We first examine the /x—>■—oo limit. Noting e = 0, we can decompose into its diagonal 

elements, which are finite, and its non-diagonal elements, which vanish: 

M^=Md + eMj (197) 

where M^i is the matrix containing the diagonal (escape) rates, and Mj is the matrix containing the 
non-diagonal (jumping) rates. The entries of Md are given by: 

L-l 

Md{C,C) = -(1 - ni)a - ^ 7x^(1 - n^+i) - 7x^/3 (198) 

Z=1 
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where is the occupancy of site i in C. At lowest order in e, those are the eigenvalues of since the 
non-diagonal part vanishes. 


Since we are looking for the highest eigenvalue of we see that there are four possible situations 
(assuming that a and fi are limited to [0,1]): 

• if a < /3, then the best configuration is empty {rii = 0 for all i’s), with an eigenvalue oi E = —a. If 
/3 < a, we have the same in reverse: the best configuration is full (rii = 1 for all j’s) and E = —/3 
(those two first cases are symmetric to one another, and we will only be considering the first one) ; 

• If a = /3 < 1, then E = —a, and we have two competing configurations: empty or full ; 

• if a = /3 = 1, then any configuration with a block of I’s followed by a block of O’s has an eigenvalue 
oi E = —1, which is the highest, and there are L -I- 1 of those. 


The phase diagram of the model in that limit thus consists of two phases, one transition line, and one 
special point (fig.-141. 



FIG. 14. Phase diagram of the open ASEP for very low current. The mean density proEles are represented in 
red the insets. The profiles in orange are the individual configurations which compose the steady state. 


However, we will be looking to perform a Legendre transform on E[^), so that this first order of the 
largest eigenvalue, which is independent of /r, will not be enough: we will have to calculate the following 
order as well. We will therefore treat the non-diagonal part of perturbatively to extract its first 
non-trivial contribution to the largest eigenvalue. 


1. Empty/full phases 


We first consider the case where a < j3, where the dominant eigenvalue of is equal to —a at leading 
order in e. We expand that eigenvalue and its corresponding eigenvector as a series in e: 

E{fM)=Eo + Y,Eke’^ ^ \Pu) = \Po)+Y.^’'\Pk) ( 199 ) 

where 

E{pi)\P^)=M^\P^). ( 200 ) 
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We already know that 


|Po) = |00 ... 00) , Eo = -a. 


( 201 ) 


The first order in s in (200) gives: 


E,\Po) + Eo\Pi) = M,|Po) + Md\Pi). 


( 202 ) 


Since the only state that can be reached from |Po) through Mj is |10 ... 00), which has no overlap with 
|Po), we get: 


\Pi) = 


1 




Af,|Po) ~ I10...00) , Pi=0 


so that the first correction to P(/i) is 0. The second order in e in (200) gives: 

E2\Po) + Eo\P2) = M,|Pi) + Md\P2). 


(203) 


(204) 


Again, the only state that can be reached from |Pi) through Mj is |010... 00), which has no overlap 
with |Po): 


lA) = ( 


1 ^2 

-M, 

Eo-Md ^ 


)Vo)--I010...00) , P2=0 


(205) 


and once more, the correction to P(/r) is 0. 

The first possible non-zero correction to P(/r) that we might get is when we find a jP^) which has an 
overlap with |Po). The shortest way to go back to |Po) through jumps is to have one particle enter the 
system (the first step being |Pi)), and then travel all the way to the other end, and exit the system. This 
can be done in L -|- 1 steps, so that = 0 for fc : 1..L, and 


El+i\Po) + Eo\Pl+i) = M,\Pl) + Md\PL+i)- (206) 

Putting those P -I- 1 first equations together, we get 

\Pl+i) = ~ |Po) + ... (207) 

and 

El+1 = (Po|Pl+i) = (208) 

1 — a 

Putting this back into P(/r), we get: 


and 


P(/r) -a + e^- - 

1 — a 


(209) 


g{j) = a-b jlog(j) - j(log(a/(l - a)) -b l). 


( 210 ) 


We may note that taking the limit oo in (149) gives the same result: 

a e'^ — 1 


Eip) = 


a -b 1 e'^ -b a 


1 1 u 

- --b -e'" = -a- 

1 -b a a 


\ — a 


( 211 ) 


which indicates that this expression for E[g) remains valid for all /i < 0. We will be more specific in 
section ED 

Finally, note that in this case, the second largest eigenvalue is —/? for e —)■ 0 (and corresponds to a 
completely full system), so that the gap between the first two eigenvalues of is finite and equal, at 
leading order, to AP = (/3 — a). 

The corresponding results for /3 < a (the ‘full’ phase) can be obtained by exchanging a with (3. 
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2. Coexistence line 


We now consider the slightly more complex case where a = j5 < 1. This time, there are two states 
with equal eigenvalues for ^ = —oo: 

|Po) = |0) = I00...00) with Eo = -a (212) 

and 

|Po) = |1) = I11...11) with Eo = -a. (213) 

As in the previous case, the first corrections to those eigenvalues are the rates with which we can go 
from these conhgurations back to themselves, but since they are degenerate, we must also consider how 
we can go from one to the other. As before, it takes the same L + 1 steps to go from |0) to itself, or 
from |1) to itself, so that the first correction to both E and E is At this stage, those states are 

still degenerate. To lift the degeneracy, we have to consider the shortest way to go from |0) to |1), or the 
opposite. This means completely hlling or emptying the system, and can be done in L{L + l)/2 steps. 
This tells us that the difference between the two highest eigenvalues is of order For 

symmetry reasons, the main eigenvector is then ^dO) + |1)), and the second one is ^dO) — |1)). 

In conclusion, we have, as in the previous case. 


A(/r) ~ -a + e^- - 

1 — a 

and 


g{j) =a + j log(j) - j(log(a/(l - a)) + l) 


but this time, the gap behaves as AE ^ e^^. 


(214) 


(215) 


To put these considerations in a more systematic format, we may use the so-called ‘resolvent formalism’ 
We present it here for two reasons: it will be useful to us in the next section, and it allows, in 
the case of a perturbation around degenerate states, to rigorously define an effective interaction matrix 
between those states. This gives us, in essence, a reduced dynamics for the system in the subset of phase 
space which contains only the dominant configurations. 

This formalism can be stated as follows: for a general matrix M with eigenvalues E^ and eigenvectors 
I Pi) and (Pi I, we may write 


i 

Jq i2TT z — M 


EiGC 


(216) 


where the sum is over the eigenvalues of M which lie inside of the contour C. Moreover, we have 


dz 


Q i27T z — M 


E E^\P.){P^\. 

EiGC 


(217) 


Going back to the matter at hand, which is with a = /3, a good way to isolate the two dominant 
eigenvectors is to consider that same contour integral, with a contour close enough to —a so that the two 
highest eigenvalues are inside it, but not any of the others. This allows us to define an effective matrix 
^eff such that: 


Meff = —a -I- 


dz 


z27r z — a — Md — eMj 


(218) 


where C is a small circle centred at 0. 

We can now expand this expression in terms of e: 


Meff =-a + 



M, 


1 


z - {Md -I- a) V z - {Md + a) 


(219) 
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which is a sum over paths of length k, with transitions given by Mj and a ‘potential’ given by (z — — 

a)~^. We see that the only terms which contribute to the integral (i.e. that give first order poles which 
yield non-zero residues) are those for which M 4 is taken at —a exactly twice, which is to say the paths 
that go through |0) or |1) twice. 

It is now fairly straightforward to find the amplitudes of Meff between |0) and |1): we only have to 
project that expression between those states, and since = —a in both of those, we only have to 
consider all the paths going from one of those states to another without going through them at any other 
point. Since we are doing a perturbative expansion in e, we only need the term with the lowest number 
of steps. Between |0) and itself, or |1) and itself, there is only one path of the lowest length, which is 
L + 1, and the amplitude for that path is y^. Between |0) and |1), or the opposite, the shortest length 
is L{L -|- l)/2. There are many suitable paths for that transition, and the total amplitude is a quantity 
X which we do not need explicitly, since that factor does not appear in the dominant term in E{^) (it 
would be of order e^^). 

All in all, we have an effective matrix given by: 


Meff = 


—a 


Xe2M 




—a 


( 220 ) 


which is easily diagonalised, and we can retrieve the results we found earlier. 


3. Equal rates point 


For the last case, where all the jumping rates are equal {a = f3 = 1), we find L -I- 1 states with an 
eigenvalue equal to —1 for fi = — 00 . Those states are given by \k) = |{l}fc{0}i_fe), i.e. configurations 
made of a block of I’s followed by a block of O’s. Those are called ‘anti-shocks’, being symmetric to the 
usual shocks which have a low density region followed by a high density one. 

Using the resolvent formalism, we find: 


{k\Meff\k) 

{k + l\Mgff\k) 

{k-l\Meff\k) 

as well as terms of the type 

{k + 2\Meff\k) 

{k-2\M,ff\k) 
{k + 3\Meff\k) 


-l + e^+\ 

(221) 

^ 5 

(222) 

L — fc-1-1 
^ ? 

(223) 


(224) 

Y^2L — 2k-\-3 

(225) 


(226) 


and so on. We can check those last terms to be of sub-leading order in and we will neglect them 

right away. 

We are left with 


L 

M^ff = -1 + +'^e'^\k){k - 1\ + e^-'^+^\k - l){k\. (227) 

k^l 

We can transform it through a matrix similarity to have all the non-diagonal coefficients be equal to 
which yields: 


L 

Meff = -l + \k){k -l\ + \k- l)(fc|. (228) 

fc=i 

This is a well known tridiagonal matrix, used for instance to model the electronic interactions in 
conjugated dienes through the Hiickel method m- It is easily diagonalised (which is left as an exercise 
to the reader). Its eigenvalues are 

Eik) = _1 + cos(fc7r/(L -H 2)) (229) 
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for fc G [1, L + 1]. The highest one is 


= -1 + cos(7r/(L + 2)) 


and the gap to the second one is 


Ai? = 2e^‘^'''^^/^^cos(7r/(T + 2)) — cos(27r/(T + 2))^ 



(230) 


(231) 


Ultimately, we find that 


E(^JL) ~ -l + 2e''/^ 


(232) 


and 


g{j) = l + 2j log(j) - 2 j. 


(233) 


Moreover, knowing that the eigenvector associated to that first eigenvalue is distributed according to a 
sine function, which is to say that the probability of |fc) is: 


P(fc) ~ sin 


nk 

L + 1 


(234) 


we find that the mean density at site n is of the form: 


Pn — 1 


+ :^sin 


L + 1 27r 


27rn 

L+1 


(235) 


Note that this probability, being given by the product of the coefficients of the right and left dominant 
eigenvectors on state |fc), corresponds to the probability of observing \k) conditioned on a low current in 
the original process M^ff, even though it is obtained from M^ff, because these matrices are similar. 


B. High current limit 


We now consider the limit where p —)■ oo. In this case, the diagonal part of is negligible, and the 
non-diagonal part is dominant. Since this non-diagonal part depends on p, we do not need to conserve 
the sub-dominant part in order to transform E{p) to g{j). To make certain upcoming calculations easier, 
we will shift our choice of pi's slightly, to 


so that we may write 


with 


Lo 
Pi 
PL 

~ {2al3e>^)T^Mj 

1 1 
^ n—1 ^ 


L + 1 

+ 

L + 1 

P 

+ 

1 

L + 1 

L + 1 

P 

_L 

1 

L + 1 


L + 1 


\og(2al5) - log(72a), 
log(2a/3), 

log(2a^) - log(v7/?), 


(236) 

(237) 

(238) 

(239) 

(240) 


where are the operators for the creation or annihilation of a particle at site n. Note that the 
dependence in p, but also a and /3, is only in a global pre-factor in M^, and hence in E{p). This tells 


49 






US that we should not expect any phase transitions in this limit (except perhaps at a = 0 or /3 = 0): the 
largest eigenvalue of will simply be that pre-factor multiplied by a constant: 


£;(Ai) (X (2a/3e'')TTT (241) 

This tells us already most of what we want to know about the behaviour of the large deviations of the 
current for extremely large positive fluctuations, but we can learn much more, as it turns out that Mj is 
exactly diagonalisable. 

We may recognise Mj to be the upper half of the Hamiltonian of an open XX spin chain m- 
Moreover, it happens to commute with its transpose. We know, from the Perron-Frobenius theorem, 
that the highest eigenvalue of that matrix is real and non-degenerate. It is therefore also the highest 
eigenvalue of its transpose, with the same eigenvectors (because they commute). This allows us to define 
H = \{Mj -I- *Mj), which has the same highest eigenvalue and the same eigenvectors as Mj. H is given 
by: 



(242) 


which is the Hamiltonian for the open XX chain with spin 1/2 and extra boundary terms S'f and iS'£ (with 
= S~^ + S~). This spin chain was studied for general boundary conditions in |174j . We will present 
here a simpler version of their calculations, only applicable to our situation but much less intricate than 
the general solution. 

The main problem in dealing with H is that it is not entirely quadratic. The first step in diagonalising 
it is to remedy this by considering two extra sites, one at 0 and one at L -I- 1, which we couple with our 
system by defining an new Hamiltonian: 


H = ^ E(5-5^+1 + S+S-^,) + ^ 

^ n—1 ^ 


L-1 


QX QX 


(243) 


Since [H, Sg] = [H, S'f+i] = 0, this modified Hamiltonian has four sectors, corresponding to the 
eigenspaces of 5'g and Since each of those has two eigenvalues 1 and —1, we can recover H by 

projecting H onto the eigenspaces of Sq and where both eigenvalues are 1: 


H — ^((Ool + (IqI)® ((0l-i-i| + (lz,+i|)-ff (|0o) + |lo))® (|0l+i) + |1l-i-i))- 


(244) 


We are now left with diagonalising H. The rest of the calculation is a rather standard approach to 
quadratic spin chains. 


We first perform a Jordan-Wigner transformation on the operators : 

c„= (n(-ir™)^- , 4= (n(-ir-)*5+, 

\m—0 / \m—0 / 

5- = (fii-ir-Acn , = ( n(-l)^™^- )cj 


(245) 

(246) 


\rn—0 


\m—0 


where rim is the number of particles on site m, with values in {0,1}. This yields fermionic operators: 


{4; — ^n,n 

Ici ci I = 0 

{Cyi, Cttj} - 0. 


(247) 

(248) 

(249) 
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The elements of H become: 


^n^n+l ~ ^h^ri+1’, 

*5^+1 “ 

= {cl - co){c\ + Cl), 
SlSl+, = (cl- cl){c[^,+cl+,), 


SO that 


L-l 


= ;^(c5 - Co)(cI + Cl) + i ^(c),c„+i + cl^^Cn) + :^(4 - cl)(c^i + CL+i). 
^ n=l ^ 


(250) 

(251) 

(252) 

(253) 


(254) 


We now perform a Bogoliubov transformation |175j on H, writing it as 

H = £o+J2^kdldk (255) 

k 

with all the ffc > 0, and the dkS to be determined. 

We want the dfc’s to be fermionic, so that [17, = £kd\, which is the equation we will now try to 

solve. We have two trivial solutions with energy 0 (called zero-modes): (cq -I- Cg) and (c^_,_^ — c^+i). For 
the other solutions, we write: 


Y{k) j 

4 = ^;r(4-co) + Xl^- 


72 




n—1 


^(fc) 

Bn^Cn + (4-1-1 + CL+i) 


(256) 


and [11, dl] = £kd\ becomes: 



1 = 2£fcfor n e [2, L - 1], 

(257) 


-Bi% - Bi% = 2f for n e [2, L - 11, 

(258) 


^(fe) ^ 2£kA[^\ 

(259) 


X(k) _^(/c) ^ 

(260) 


^(fc) ^ 2£kX^'^\ 

(261) 


y W + 47^ = 

(262) 


_y(fc) = 2£kB^^\ 

(263) 


Al'’ - Bl'’ = 2£kY‘'’^\ 

(264) 

All those equations can 

be written in a more compact form by defining: 



/t(^) _ yik) 

^L+1 ~ ^ ’ 

(265) 


A ( 1 

^L+l+n ~ V ^L+l-n^ 

(266) 


II 

So 

(267) 

for which they become: 

Ai^l, + Ai% = 2£kAi^^ for n G [1, 2L1, 

(268) 


life) 1 ( Yk) __ nc Yk) 

^2L V ^0 ~ ^^k^2L+l^ 

(269) 



(270) 


These are the same equations which we would have found for a periodic XX spin chain with 2L 2 
sites, the only (but important) difference being that d\ mixes c^’s and cl’s, so that the total spin is not 
conserved. 
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We look for plane wave solutions of the form An = r", with 28 = r + ^. This automatically solves 
eq.(268|. The other two equations become: 


+ (- 1 )^ = 

1 
r 


-Si 

‘2L-\-\. I _ , 


+ r = r + 


and both simplify into 


2L+2 _ 


= (-ir 


(271) 

(272) 

(273) 


iir{L-2k+2) 

We have 2L + 2 solutions to this equation, given by r = Wfc = e 2 i ,+2 for k G |1,2L + 2], so that 
Ait'’ = ujk, and the energies are given by: 


gfc=cos( ) for fc€ [l,2L + 2]. 


2L + 2 


2L + 2 


We can then write the d^’s as: 


^ 72 LT 2 ( ^ ^ ~ (-a;fc)-”c„ + -^(c|^+r + cl+i) 


and the inverse relations as: 


1 


2L+2 


E^: 


4 


V2^^° V2L + 2 ^ 


c'l' = 


1 


2L+2 


V2LT2 f 


E ”4. 


2L+2 


= vlFfI S 


^(cl+ 1+^f-+i) ^2L + 2 


fc=l 

2L+2 

E 

fc=i 


CJ 




(274) 

(275) 

(276) 

(277) 

(278) 

(279) 


Note that the cij(,’s are fermions, but, because there are 2L + 2 of them, and only L + 2 of the cj(’s, 
they are not all independent: we have uJL+i+k = —^k, so that d^_|_r_|_j, = —dfc- 

We now need to determine the constant term 8o in ff. Considering only the scalar terms in ^ E^d^dk, 


we find: 


L+i I fI ^ 1 \ 

E 2 L + 2 “ 4 ) + E 4cn + Cnci + 2(4+1 + CL+l)(cL+l + c[_^r) 1 

k^l \ n^l / 


L+1 




(280) 


k=l 


Since there is no scalar part in H, we must therefore have: 

^ Z/+1 

^0 = ~ 2 E! ■ 

fe=l 


(281) 


The highest eigenvalue can then be obtained by considering the state with all the energy levels occupied: 


( 282 ) 
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The corresponding eigenstate is defined by 


L+l 

iV')=n 4 if^) 


k=l 


( 283 ) 


which is such that d|,|'0) = 0 for A: G |1,L + 1]. The vector \^1) is arbitrary (provided that it is not in 
the kernel of any of the d^’s that we apply to it). 

Remembering the global factor which we took out of at the beginning of this section, we finally 
get: 


TT 


(284) 


and 


aU) ~ L[j log(j) - j{l - log(7r))) 


(285) 


which is proportional to L, consistently with eq.(189) which was obtained for positive fluctuations of the 
current in the MC phase. 


We now have what we wanted, but we can have a look at the eigenvector (2831 as well. 


The first and easiest calculation that we can do here is that of the two-point correlations in 
connected correlation between the occupancies of sites n and m is given by: 


Cnm, — 


cJ,Cnc|„Cm| 


- i'tPKCnWii’lclaCmW 


= -{i’\4i4nW{i’\CnCm\'lp) + (V'l | V”) (V” |Cn4 


The 


(286) 


(using Wick’s theorem). We find that: 


(V’l4cLlV’) 

{lp\CnCm\4’) 


L+l 


—y 

, 4- 9 ^ 


2L + 2 ’ 

1 L-\-\ 

-V(-a;fc)”’ 

2L-t2^^ ’ 

k=l 

^ L+l 


2L + 2 
1 

'2L-1-2 




L+l 


yi -^ k ) 


n m 
^k 5 


k=l 


(287) 

(288) 

(289) 

(290) 


If n and m have same parity, each of those terms sum to 0 (unless n = m in the first two sums, but here 
we consider two different sites). If not, we get: 


so that 




\LnLn 


2^ gi7r(L)(m—n)/(2L-|-2) 

L -I- 1 I — ’ 

1 piTr(L){—m—n)/(2L+2) 

L-f \ — Q^'^(m+n)/{L+l) ■ 

1 rJ,'J^(L)(m+n) / (2L+2) 

_—(- 1 )" ^ 

L-fl' ' I — Q-iTr(m+n)/(L+l)' 


c ^ 


1 


4(L-hI)2si+ 

7r(m+n) 

1 (2L-e2) ) 

4(L-Hl)2si+ 

7r(m—n) 

( {2L+2) ) 



(291) 

(292) 

(293) 


( 294 ) 
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The correlations are therefore exactly 0 for sites which are an even number of bonds apart (as is the 
case for a half-filled periodic chain |104) b and behave as 


C„ 


1 

Tr‘^{m — nY 


(295) 


otherwise, if the two sites are far away enough from the boundaries. Note that those correlations do not 
vanish with the size of the system, in contrast with the steady state of the ASEP at /r = 0, where they 
behave as L~^ in the maximal current phase and vanished exponentially in the high and low density 
phases [TUg] , 


We can now examine the (un-normalised) probability of any given conhguration. This can be expressed 


as 


P{{hn,Pn}) = 


N 




L+l 


+ cl+i)' 


n =0 


L 

n 




w 


(296) 


which is to say that we select only the configuration that has holes at positions {/i„} and particles at 

positions {pn} in nnd project it onto its Hermitian conjugate. Note that the term .. . makes 

no difference (since it is equal to a constant factor ^), but is there to h ave e ffectively L + 1 sites instead 
of L, which will be useful shortly. From the expression of this term in (279), we see that it corresponds 
to having a hole (or, in fact, a particle) at site L-|- 1. Note also that the terms in (296) can be reordered 


as long as any pair {c„, c])} is kept in the same order, so that we can regroup all the c’s from the first 
product to the left, for instance. 

We can now use Wick’s theorem m on this expression, and write it as the Pfafhan of an anti¬ 
symmetric matrix A whose upper triangle consists of all the mean values {i/jlablip) where a and b are two 
terms from the product in (296), taken in the same order. We can write it as a block matrix: 


A = 


Ai 


-{Y\ch„ 

-{Y\ch„ 

-{Y\ch„ 


where Ai\n,m = {'Y\ch„Chrr.\'>P) H n < m 
with c)j , A 3 with Cp^ and A 4 with . 



4) 

(V’|c/«„4^ 

4) 

{Y\ch 

n ^Pr 


^2 



4) 

414 


J4 


J4 

A 3 


414 

-i ^Pn 


-(v-kLcp 




A 4 


and -Ai — 


4) 

if n > 

m. 

The 


(297) 


Looking at the expression given in (287) to (290), we see that —{ip\cmC-^ 
and -{iplcmcilip) = 


l4cL 


T, 


- S„ 


= {Y\CnCm\Y), -(V'|cL4lV’) = 

We also note that all those block matrices can be 


factorised in a simple way: if we define 

44 = 


-hn 


^n,k = W.. 

Y+ — 


then A can be rewritten as 


A = 


x-{x+Y x-{x-Y A-(r-)t 

A+(A+)t-l,, A+(A-)t A+(r-)t 

y+{x+y Y+ix-y Y+iY-y 

y-(x+)t r-(A-)t Y-{Y-Y-ip 


x-{y+Y' 

x+(y+)t 

y+(y+)t 

r-(y+)t 


(298) 

(299) 


(300) 


where Ih and Ip are identity matrices whose respective sizes are the numbers of holes and particles (one 
of which is occupying site L -I- 1, so that the sum of the two numbers is L -I- 1). 

In order to calculate P({/i„,p„}) = Pf [A], we need to first consider its square P({/i„,p„})^ = Det[A]. 
After exchanging a few lines and columns in that determinant, we can write it in a factorised form: 


Det [a] = Det 


y- 

- 1 l+i 

y+ 

0 


(x+)t (y-)t 

(A-)t (r+)t- 

1l+i 

0 


(301) 
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Each block in this last expression is of a square matrix of size L + 1 (which is where the fact that we 
included site L + 1 becomes useful). That determinant then reduces to: 


Pi{K,Pn})^ =Bet[^;]Bet[iX-)^ (r+)t] = |Det[^;]|^ 


After a few final simplifications, we find 


Pi{hn,Pn}) = Det 




(302) 

(303) 


where uj = e^+i. We recognise this to be Vandermonde determinant which gives, for any configuration: 


P{{nn})= n 


[sin(rm - r„)] [sin(r 

rin^rim 


(304) 


where n„ is the occupancy of site n, and r„ = mTj{2L + 2). Note that all these probabilities are still 
un-normalised. 

This distribution is exactly that of a Dyson-Gaudin gas HZZl, which is a discrete version of the Coulomb 
gas, on a periodic lattice of size 2L + 2, with two defect sites (at 0 and L + 1) that have no occupancy, 
and a reflection anti-symmetry between one side of the system and the other (fig.-15). The first (upper) 
part of the gas is given by the configuration we are considering, and the second (lower) is deduced by 
anti-symmetry. The interaction potential between two particles at positions r„ and is then given by: 


V{rn,rm) = -log(sin(rm - r„)). 


(305) 


It was shown in |104j that the large current limit of the steady state of the periodic ASEP of size L 
also converges to a Dyson-Gaudin gas (the standard periodic case, without defects or symmetry). 



FIG. 15. Dyson-Gaudin gas equivalent for the configuration (110101000110111) for the open ASEP conditioned 
on a large current. The lower part of the system is deduced from the upper part by an axial anti-symmetry. 


We should note that the trick consisting in taking the sum of Mj and its transpose to reconstruct 
an XX spin chain is not in fact necessary |178) . All the calculations we saw can be performed, in a 
slightly different way, on Mj directly, which has the added advantage that the imaginary part of the 
other eigenvalues is not lost. 


C. Generic dynamical phase transition 

As we stated in the introduction to the present section, none of the methods we used to analyse the 
/i —>■ ±oo limits appear to rely on the ASEP being integrable. It could be, however, that they do require 
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it in a subtle way. This is, in fact, not the case, as it turns out that those two limits can be treated in 
the exact same way, and yield very similar results, for a relatively broad class of models which are in 
general not integrable |179j . 

Let us consider a generalisation of the open TASEP, with site-dependent jump rates pi, as well as 
an arbitrary (but finite) extra interaction potential V(C) which may depend on the whole configuration 
C, so that the rate for a particle to jump from site i, turning configuration C into C', is of the form 
The methods we used for the TASEP can be applied in precisely the same way to that 

new model. 

In the low current limit, we found for the TASEP that the first order in perturbation of the largest 
eigenvalue did not depend on the size of the system, in every case. This comes from the fact that, whether 
the dominant eigenstate at order 0 is degenerate or not, the fastest way for the system to leave one of 
these states and come back to it is to have one particle jump through the whole system once, which 
corresponds to one quantum of current. This remains true for the generalised model we are considering 
now, so that the large deviation function of the current does not depend on the size of the system for low 
currents. 

In the large current limit, the similarity is even stronger, as was noticed in m- The non-diagonal part 
of the deformed Markov matrix which we are keeping in that limit is entirely equivalent to that of the 
simple TASEP, up to a matrix similarity, so that its whole spectrum is the same (up to a multiplicative 
constant, in fact). The eigenvectors are also related to those of the TASEP, and the product of the 
right and left dominant eigenvectors, which give the stationary density conditioned on a large current, is 
identical. More importantly, the large deviation function of the current is linear in the size of the system 
for large currents. 

One will be able to find a detailed account of these statements in mi- 

In the case of the open TASEP, these two limits are compatible with the results obtained for small 
fluctuations of the current around the MC phase: the large deviation function is independent of the 
size of the system for negative fluctuations, and proportional to it for positive fluctuations, so that a 
dynamical phase transition takes place at the average current. In the next section, we will see that this 
comes from a transition in the appropriate coarse-graining scheme in the large size limit. For currents 
which are below |, which is the average current in the MC phase, the system can be describes through 
a hydrodynamic (i.e. mean field) description, where the large size behaviour of the system depends only 
on the local density. In that scenario, as we shall see shortly, the optimal way to produce a fluctuation 
of the current is to have a localised fluctuation of the density. The localised nature of that fluctuation 
means that its cost, in terms of probability, will not depend on the size of the system. 

However, as we saw in section |III B[ the largest current which can be produced in the mean-field 
approach is In order to have the current fluctuate above that value, one has to introduce correlations, 
spread out through the entire system, which can be seen in the large current limit. Since these correlations 
are needed everywhere in the system, the probability cost is extensive, hence the factor L in the large 
deviation function of the current. Moreover, because of those correlations, the local density fluctuates 
much less that in the hydrodynamic case, so that these states are sometimes called ‘hyperuniform’ [93]. 

One would hope that these observations remain true in the case of generalised rates, and that the 
behaviour of the large deviation function of the current in these two limits is a good indication of the 
presence of a dynamical phase transition somewhere in between. It is not known, for the moment, whether 
there are any universal features to that transition. It would be particularly interesting to examine the 
case where the jumping rates are disordered. 
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VI. HYDRODYNAMIC DESCRIPTION AND DYNAMICAL PHASE DIAGRAM 


In the previous sections, we obtained, by various exact calculation methods, the behaviour of the 
large deviation function of the current of the open TASEP in three different limits: around the average 
current, and for very high and very low currents. In this final section, we will see how the gaps between 
those limits can be bridged through a hydrodynamic description of the system, based on the macroscopic 
fluctuation theory (MET, [73]). Although this method is much less rigorous than the exact calculations 
performed earlier (because it relies on a poorly controlled exchange of limits), we will see that it gives 
much more information on the fluctuations of the currents, in the cases where it applies, and gives us 
access to the complete dynamical phase diagram of the open ASEP. Five phases will be obtained: a low 
and a high density phase, which are the continuation of those same phases for the average current, and 
which meet the empty and full phases in the low current limit ; a shock phase, which stems from the 
shock line, for positive fluctuations of the current ; an anti-shock phase, which is obtained by decreasing 
the current from the maximal current phase of the average current ; and a maximal current phase, which 
is obtained for j > \. In the first four of these phases, the large deviation function of the current is 
independent of the size of the system, and the results obtained with this method agree with what we 
obtained earlier in the appropriate limits. In the last phase, the hydrodynamic description breaks down, 
as the system goes through the dynamical phase transition which we discussed at the end of the previous 
section. 

We will start by presenting the principle of the method, which was first used in [81Ij on the open 
ASEP, but works equally well for more general one-dimensional driven lattice gases. We will then use it 
to describe all the dynamical phases of the ASEP, obtaining the large deviation function of the current, 
as well as the form of the typical profiles conditioned on that current, save for the maximal current phase 
where the method breaks down. This will allow us to draw the dynamical phase diagram of the system. 
Finally, we will comment on evidence showing that this method yields not only the dominant state and 
eigenvalue of the deformed Markov matrix, but also a large family of higher eigenstates. 


A. Macroscopic Fluctuation Theory: from WASEP to ASEP 


The method which we will be using in this section, and which can be found in [S3], is based on the 
macroscopic fluctuation theory m, or MET for short. The principle is the following: for a noisy Langevin 
equation with a Gaussian noise, such as a noisy Burgers equation, the probability of a certain history is 
given by the Gaussian weight of the noise which generates it. For instance, for a driven one-dimensional 
interacting particle system, with a conductivity cr[p) and a diffusion coefficient D{p)^ subject to a field 
F, the local density verifies the equation 

dtp = -Vj with j = Fa{p) - D{p)Vp + (306) 


where ^ is a Gaussian noise with mean 0 and variance 1, and it is implied that p and j depend on x and 
t. The probability of a certain history p(x,ty\ is thus given by 


^og(p[{j{x,t),p{x,t)}]'^ dr 


[j - Fa{p) +D(p)^p\‘ 

2a{p) 


•dx 


(307) 


with the constraint that AtP = This gives explicitly the joint large deviation function of the current 

and density. In order to obtain that of the current alone, one then has to minimise that quadratic action 
with respect to the density (as stated by the contraction principle ; c.f. section IIA). This minimisation 
will give us not only the large deviation function of the current, but also the optimal density profile to 
produce that current, which is the value of p realising the minimum. 

However, this cannot be directly applied to the ASEP: looking back at to eq.(63l, which is the de¬ 


terministic part of the Langevin equation, we see that D = vanishes at large sizes, or, equivalently 
(through a rescaling), that the field F ^ L(1 — q) diverges, leaving us with an inviscid Burgers equation. 

We will use this approach nonetheless, by starting from the weakly asymmetric simple exclusion process 
(or WASEP), for which F is finite (which is to say (1 — 9 ) ~ L~^), or in fact the general system described 
by eq.(306). We will then minimise the MET action with respect to the density, after which we will take 


F to be of order L again. 
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We therefore start with the joint large deviation function of the current and density, given by 


9{j,p) 


lim - 

t —^-j-CXD t/ 



[j - Fajp) + J(p)Vp] 
2 cr(p) 


2 

-dx 


( 308 ) 


with dtp = —Vj. We intend to contract this to the large deviation function of only the space integrated 
current. We will assume that the best density to produce a given constant average current, in the long 
time limit, is time-independent. The constraint thus becomes dtp = 0 = —Vj, which is to say that j 
does not depend on a; or t any more. The time integral can then be taken out of the large deviation 
function, and we get 


9{j,p) 


/•i [j - Fa{p) + D{p)Vp]\ 

Jo Mp) 


(309) 


which is the action we will minimise over p. 

The first step in that direction is to expand the square 


9{j,p) 


'■i [j-Fa{p)]+[D{p)Wpy 
3 2cr(p) 


•dx + 


[j - Fajp)] 

^(P) 


D[p)\Jp dx 


(310) 


and notice that the cross product between the gradient and the rest produces a constant term 


] ^[P) jp^ cr(p) 


(311) 


which does not need to be minimised. We are left with having to cancel the functional derivative of the 
first part alone. 

Let us write 


A{p) 


[j - Fa{p)]^ 
2 cr(p) 


and B{p) 


[Djpf 

2 cr(p) 


(312) 


We want to minimise f^dx A(p) + B(p)(Vp)^. The minimising profile satisfies the Euler-Lagrange 
equation: 

A'{p)-B\p){Wpf-2B{p)Ap = 0 (313) 


which, multiplied by Vp, gives 

V[^(p)-i?(p)(Vp)2] =0. (314) 

The minimising profile is therefore such that 

A{p)-B{p){Wpf=K (315) 

where K is an integration constant. That constant can be found to be 0 for E —>■ oo [SS]. Replacing A{p) 
and B{p) by their expressions, and taking a square root, we obtain the equation defining the optimal 
profile: 


J - Fa{p) ± D{p)Vp = 0. 


(316) 


This equation is the same as the mean field equation, except for the sign of the gradient. To construct 


this profile, for a given current, one may therefore use the same construction as we saw in section IIIB 
but allow the gradient to have either sign. We will come back to this in the case of the ASEP. 


We also obtain, from this minimisation, the large deviation of the current, by injecting eq.(316) into 


p(j, p). The portions of the optimal profile where Vp has the same sign as in the mean field equation do 
not contribute to the integral, as they make the integrand vanish, so that: 


2 cr(p) 


2ct(p) 


<P) 
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where the integral is over all the portions of space where j — Fa{p) = D{p)Vp. We end up with a 
combination of the values of the primitive of (j — Fa{p))/a{p) at the points delimiting those portions of 
space. 

We now consider the case of the TASEP, where F = L, a{p) = p{l — p) and D = ^, and we rescale 
time by a factor L, so that the optimal prohle equation for a current j is given by 

j = p{l -p)T ^Vp. (318) 

For a given current j, and two boundary conditions Pa and p^, we can repeat the construction of density 
profiles we saw in section [IlIB[ without the constraint on the sign of the variations of p. Note that this 
constraint enforced that for given boundary conditions, only one value of the current was possible. In the 
present case, the current is a free variable. The allowed building blocks for the profiles are represented 
on the portions in red are the only ones contributing to g{j). 



FIG. 16. Various optimal profiles for a fixed r. The parts in green satisfy the mean field equation, and do not 
contribute to g{j), whereas the portions in red do. 


The large deviation function of the current for that profile is then given by an integral over the red 
portions of the profile: 


gij) = 


P(1 -P) ^ 


P 


I- p 


1 Pi 


(319) 


where p^ and pf are the densities at the boundary of each red section. Note that in every case, one of 
the boundaries is r or 1 — r. We shall write each of these terms as 


/(j;rl,r2)=^ dp = j log( 


' 1 — n r 2 ' 
ri 1 - r2. 


+ fi - r2 


(320) 


where j has to be equal to either r’i(l — ri) orr 2 (l — r 2 ). This is the same as the function E’'®® from [53] . 

Through this procedure, there are many profiles we can build for a given set {j,Pa,Pb\, since any 
number of shocks and anti-shocks can be added between r and 1 — r. The true optimal profile is the 
one which minimises the large deviation function, which is to say the one with the fewest red portions. 
If at least one of r and 1 — r is between pa and pf,, then the optimal profile is monotonic. Otherwise, 
two candidates must be compared: one where p goes from pa to r, and then to pb, and one where p goes 
from ptj to 1 — r, and then to pb- 

The other extremising profiles are conjectured to correspond to excited states of the system, condi¬ 


tioned on the average current. We will say more about these in section VIC 
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In the next section, we will list all the possible configurations for j = r(l — r), pa and ph that lead 
to different forms of g{j) (i.e. with different combinations of the function /), and determine, in each 
of those phases, the expressions of ^(j), E{p), j{p,), and the boundaries of the phase. We will also 
compare the asymptotic behaviours of those results with everything we found in the previous sections, 
to confirm their validity. We will then summarise all we know about the maximal current phase, which 
is not accessible by this method (and which is defined by j > \). Finally, we will put all this together in 
order to draw the phase diagram of the open ASEP with respect to pa, pi, and p,. 

In all cases, we will be noting r the density for which j = r(l — r) which is below i. Also note that we 
will do all the calculations for the TASEP, knowing that the same for the ASEP can be obtained merely 
by multiplying A(p) by {1 — q). We also recall that we have defined two other boundary parameters 
a = and b = which we will use in certain formulae to make them more compact. Finally, we 

will, for the same reason, sometimes use, instead of /i, the variable u defined by: 


1 

1 +e'^’ 


(321) 


The non-perturbed case is given by it = ^, it = 0 corresponds to an infinite current, and it = 1 to zero 
current. 


B. Dynamical phase diagram of the open ASEP 

Each phase is described in turn. 


1. High/low density phases 


We start with the low density phase, from which we can deduce the high density phase through 
Pa^l- Pb- 

This phase is defined by Pa < 1 ~ Pa < 1 — i" and pf, < 1 — r. The optimal profile is almost always on 
p = r, with possible boundary layers at both ends, with only the one at the left boundary contributing 
to g{j) (fig.-^l. 



FIG. 17. Optimal profiles for a fixed pc in the low density phase. Only the portion in red contribute to g{j)- 


The large deviation function of the current is in this case: 


gU) = fi.r,Pa,r) = jdog( 


1 - Pa r 
Pn. 1 - r 


) + Pa - -r, 


(322) 


which agrees with what we found in eq.(150l from exact calculations. This was first obtained in 
The generating function of the cumulants of the current is 


E{p) = 


a — 1 


1 


a + le'^ + a e^^ + a l + o’ 


(323) 
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which was also obtained in |124) through a numerical approach to the Bethe equations, and the current 
is, in terms of 


3 ip) = 


a e^ 


(e^ + a)2 

The boundaries of the phase are given by: 

Pa <1- Pb 


= Pa(l - Pa) 


u{l — u) 


{pa+U- 2upa)'^ ' 


U > 


U > 


pI 


l-2pa+ 2pl 
PaPb 


with Pa> 


1 ,0 - with Pb > 

1 Ph Pa “t” 2paPb 2 


U> Pa 


1 1 

With Pa < 2 > Pb< 2’ 


(324) 

(325) 

(326) 

(327) 

(328) 


where this last condition corresponds to j < 


j, which is the boundary with the MC phase. 


According to this, the LD 
with what we found for p —>■ 
We may also note that, on 
we find: 


phase goes all the way up to u = 1. This expression of E{p) is consistent 
—oo, i.e. u —>■ 1, i.e. eq.(210). 

the line Pa = ht which corresponds to the LD-MC transition line for p = 0, 


E{p) 


le'" - 1 
2et^ + l 


(329) 


which is consistent with the expression found in [97] for the half-filled periodic TASEP (we recall that 
the open system with Pa = \ and pb < 1/2 is equivalent to a half-filled periodic system of twice the size). 
The p —^ 0“ limit gives: 


E{p) 


p _p^ 

4 48 


which is the same as eq.(179|. 


(330) 


2. Shock phase 


We consider the case where Pa < r and pb > 1 — r. Here, there is a number of optimal profiles of 
order L. Each of them has a boundary layer around each boundary, both of them contributing to <?(/), 
and two constant regions, where p = r near the left boundary and p = 1 — r near the right boundary, 
separated by a shock that can be placed anywhere in the system (fig.-18). 



0 L 

FIG. 18. A few optimal profiles for a fixed Pc in the shock phase. Only the portions in red contribute to g{j). 


The large deviation function of the current is given by: 

'il-pa)pb r'^ 


aij) = fir^Pa, r) + /(/; 1 - r, pb) = j log 


Pa(l - Pb) (1 - rf 


+ Pa - Pb + 1 - 2r. 


(331) 
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The generating function of the cumulants of the current is 

2 e '"/2 I 


E{p) = 


e'^/^ + ^/ab 1 + a 1 + & 


and the current is: 


jXm) = 




Pa)pb{l-u) {l-pa)pb , /(1-u) 


(eM/2 _|_ y p^(i _ pj) 

The boundaries of the phase are given by: 


Pa(l - Pb) 


( 332 ) 


(333) 


u < 


u < 


u > 


PaPb 

with 

Pa 

1 

<2 ’ 

Pb>l, 

(334) 

1 Pb Pa ^PaPb 

(1 - Pa)(l - Pb) 

with 

Pa 

A 

to 1 

Pb>l, 

(335) 

1 Pb Pa H" ‘^PaPb 

Pa(l - Pb) 

Pb Pa ‘^PaPb 

with 

Pa 

1 

^ 2 

, Pb>l, 

(336) 


where the last condition corresponds to j < We may note that the volume which is defined by these 
boundaries is symmetric under any permutation of pa, 1 — Pb and u. 

The shock phase concerns one of the asymptotic results we found before. For p —>■ O'*', which imposes 
Po = 1 - Pb, we find: 


E{p) 


a a{a — 1) 

(l + a)2^^ 4(a+l)3' 


(337) 


which is what we found in eq.(159). 


3. Anti-shock phase 


The last phase we can access through the MFT is for pa > (1 — r) and pb < r. In this case, there also 
is a number of optimal profiles of order L: the first go down from p^ to (1 — r), then down from (1 — r) 
to r through an anti-shock that can be placed anywhere, and that contributes to g{j), and finally down 
from r to pb (fig.-191. 
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Pa 
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0 


1 



Pb 
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L 


FIG. 19. A few optimal profiles for a fixed pc in the Anti-shock phase. Only the portions in red contribute to 

gU)- 


The large deviation function of the current is given by: 

gU) = /(j; 1 - r, r) = 2j -h 1 - 2r. (338) 

The generating function of the cumulants of the current is 

2et‘/2 

= 3^/2 1 = tanh(p/4) (339) 
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and the current is: 




1 — tanh(/r/4) 
4 


—2{u — V?) + y/u — v? 
1 — 4(tt — u^) 


The boundaries of the phase are given by: 


u < 


u < 


u> t: 


1 - 2p„ + 2pl 

1 ~ ‘2‘Pb + ‘2'P^ 

1 

2 ’ 


with Pa> ^ 
with Pb 


) Pa > 1 Pfc) 
; Pa ^ 1 Pb; 


( 340 ) 


(341) 

(342) 

(343) 


where the last condition corresponds to j < j. 

We note that this phase corresponds to one of the examples that can be found in [S3] . The expression 
for E{p) also comes up as a side note in [124j . 


The limit /r —>■ 0 gives: 


"" 4 192 


which is consistent with eq.(193|. The limit p, —?► —oo, which implies p^ = 1 — p;, = 1, gives: 

Eip) - -1 + 2e'"/^ 


(344) 


(345) 


which is the same as equation (232), and this is the last asymptotic limit that we had to check. 


4- Maximal current phase 

There is one phase left for us to examine, to a much lesser extent than all the the others because the 
MFT breaks down in this case: the maximal current phase. Once we take out the phases we have already 
considered, we are left with a volume, in the three-dimensional phase space with variables pa, Pb and u, 
defined by: 


1 

"<2 
U < Pa 


U < 1- Pb 

Pa(l - Pb) 


with Pa> 
with Pa > ^ 


with Pa < 
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P.< 2, 


Pb> 2. 


Pb < 


1 

2' 

1 

2 ' 

1 

2 ’ 


u < 


with Pa < 

Ph 4” Pa ^PaPb 2 


P6> 2- 


We know that, asymptotically: 


P(j) (J - J)'/" 


32v^L 
57 r(l — (7)3/^ 


(346) 

(347) 

(348) 

(349) 

(350) 


for p —)• 0+ with Pa > \ and pb < ^ (i.e. right next to the MC phase for the steady state), which we 
found in eq.(175), and that: 


g{j) - Lj log(j) - Lj{l - log(7r)) 


(351) 


for p — 00 , as we saw in eq.(285). The best estimate of g{j) from the MFT is obtained for p = 1 in the 
whole system, which yields g{j) ^ L{j — J)^, showing that this hydrodynamic description does indeed 
break down in the MC phase. 
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FIG. 20. Phase diagram of the open ASEP in the s-ensemble. Top: diagrams at fixed u. Centre: complete 
diagram with phase boundaries and exploded view. Bottom: diagrams at fixed pa- 
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Since this last result is valid independently of the boundary parameters, we know that all the plane 
u = 0 belongs to the same phase. We have, however, no way to be certain that all of the volume we have 
described above is just one phase. That being said, we know that the entire region corresponds to a mean 
current higher than ^. There is no way for the system to produce such a current through a hydrodynamic 
profile, for which the maximal possible current is ^ if p = |, so that, in order to increase the current, 
the system must produce correlations, which is why the MFT breaks down. Those correlations must be 
negative for neighbouring sites (if the particles are next to holes, they will jump more easil y an d produce 
more current), which is consistent with what we found in the large current limit in eq.(295). What’s 
more, those correlations must be created everywhere in the system, because a single non-correlated zone 
would cause a blockage and bring the current back down to ^. We can also argue that the mean density 
should be around because it is easier to get to a large current starting from ^ than from anything 
lower. In all these remarks, the boundaries play little part: independently of them, the system must 
be around p = f, and anti-correlated at every point. We therefore don’t expect any sub-phases in this 
region. 


5. Phase diagram 


Now that we have considered all the possible combinations of p^, pb and u, we can draw the phase 
diagram of the ASEP in the s-ensemble (fig.-201. Each phase is represented using a different colour: 


blue for the low density phase, green for the high density phase, orange for the shock phase, purple for 
the anti-shock phase, and red/pink for the maximal current phase. The full diagram can be seen in the 
centre of the figure, with black lines marking the corners of the phases, and an exploded view of the 
LD, MC, shock (S) and anti-shock (AS) phases is also shown, with coloured lines representing slices for 
regularly spaced values of u. The HD phase can be deduced from the LD phase through the symmetry 
Po O 1 — pt,. The top and bottom parts of the figure contain slices of the diagram for specific values of 
u (top) and Pa (bottom), with a few iso-current lines drawn in all phases except the MC phase. Note 
that those iso-current lines do not represent evenly spaced values of the current (j varies, in fact, more 
slowly as one approaches the MC phase). 


C. Hydrodynamic excited states 


As we noted earlier, the extremisation of the MFT action yields many density profiles, which are local 
minima of that action, and of which we only considered the optimal ones to build the dynamical phase 
diagram of the current. One may naturally wonder if the other extremising profiles have any physical 
significance. We conjecture that they do in fact correspond to the low excitations (i.e. relaxation modes) 
of the model. In particular, profiles which correspond to a vanishing p, i.e. ^g{j) = 0, if they exist, are 
relaxation modes of the unperturbed dynamics, and the corresponding values of the Legendre transform 
of g{j) are the eigenvalues of these relaxation modes, which is to say the inverses of the relaxation times. 
In this section, we give some evidence towards that conjecture. More details on this conjecture will 
appear in m- 


Consider the system in its low density phase, with some fixed values of j = r(l — r) and pa-, and 
pi, < I — r. Looking for the second best profile obtained from the aforementioned extremisation, we 
notice that it depends on the position of pb with respect to r: if pb < r, it is obtained by adding a shock 
and an anti-shock to the optimal profile, whereas if pb > r, one may go directly down to pb after the 
shock, so that the profile with an anti-shock becomes the third best (fig.-211. 


In the case where pb < r, the large deviation function of the current for the second best profile is given 


by 


gU) = ^ (r^) ^+i + Pa-3?'- 


The deformation parameter p vanishes for 


1 ^ ^ 1 / 31-1 

Pa 


(352) 


(353) 
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FIG. 21. Second best extremising profiles in the LD phase, for pt, smaller or larger than r. The profile on the 
left still exists for ph > r, but becomes the third best. 


and we get a value for the Legendre transform of g{j) equal to 

E* = —1 — Pa + ?>Pc- 


(354) 


If Ph > r, that profile still exists, but another one appears, which is more probable (fig.-21 right). The 
large deviation function of the current for that profile is 


gU) =j^og(^—Pa Pb , Pc ^2\ 2 p^ 

V Pa I - Pb^l- Pc J 


Pa 1 Pb 1 , 

The deformation parameter p vanishes for 

'1-Pa Pb \i/2l“i 


1 + 




Pa 1 Pb 

and we get a value for the Legendre transform of g{j) equal to 

E* = — Pa + Pb -\- 2pc. 


(355) 


(356) 


(357) 


1-pb 



FIG. 22. Phase diagram of the first excited state of the open ASEP ; an extra transition appears in the LD and 
HD phases, corresponding to a change of behaviour of the associated density profile. 
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The line on which the second profile appears is that for which 


Pb 




which is to say, using the alternative boundary parameters a and b: 


( 358 ) 


ab-^ = 1. 


(359) 


This line is therefore a transition line in the phase diagram of the second best density profile of the open 
ASEP (hg.- 221. The same calculations can be done for the HD phase by exchanging a and b. 

The value obtained for E*, as well as the location of this transition line, agree perfectly with results 
obtained in [35] for the eigenvalue of the first excited state of the open ASEP, through a numerical 
resolution of the Bethe equations, and which were later verified in |123j . This leads us to conjecture that 
all the non-optimal extremal profiles obtained from the MET are excited states of the (biased, unless 
/i = 0) open ASEP, with the value of the Legendre transform of the large deviation function giving the 
corresponding eigenvalue. One should note that each shock and anti-shock in the profile adds an order L 
to the degeneracy of that eigenvalue (because it does not depend on their position), and that finite size 
corrections will lift that degeneracy, in which case the true eigenstates will be specific superpositions of 
these profiles. More evidence for that statement will appear in [181] . 
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VII. CONCLUSION AND OUTLOOK 


We have seen, in this review, how to obtain the large deviation function of the average current of parti¬ 
cles in the steady state of one-dimensional bulk-driven lattice gases, and in particular of the asymmetric 
simple exclusion process with open boundaries. 

We hrst reviewed the mathematical tools necessary to pose and treat the problem at hand. We defined 
large deviation functions, which are a natural generalisation of free energies, and saw that they relate to 
generating functions of cumulants through a Legendre transform. We then considered the special case of 
time-additive observables in time-continuous Markov processes, and saw how the problem of obtaining 
the generating function of those observables in the long time limit reduces to that of computing the 
largest eigenvalue of a deformed Markov matrix. 

We then introduced the reader to the open ASEP, starting with a rapid overview of the existing 
literature related to the model and its variants. We also presented two important results pertaining to 
its steady state: first, the phase diagram of the average current, as well as the typical density profiles, 
using a mean-field approach ; then, the famous matrix Ansatz, giving the exact probability distribution 
of the steady state using a matrix product formulation, and which allows to obtain that same phase 
diagram through an exact calculation. 

Being interested in the fluctuations of the current rather than its average, we posed the problem of 
obtaining its large deviation function, which is equivalent to that of calculating the largest eigenvalue 
of a deformed Markov matrix. That matrix being integrable, we saw how to use the coordinate Bethe 
Ansatz (in the periodic case) and the Q-operator method (in the open case) in order to obtain an exact 
expression of the generating function of the cumulants of the current, written as an implicit pair of series 
in an intermediate parameter. Treating those series in the large size limit, we obtained the asymptotic 
behaviour of the large deviation function of the current for small fluctuations, in each of the phases of 
the average current, and found a few dynamical phase transitions. We noticed in particular that for 
fluctuations of the current smaller than j , which is the maximal current obtained from a hydrodynamic 
(mean field) description of the system, the large deviation function does not depend on the size of the 
system, whereas if the current is larger than j, it is proportional to the size of the system. 

Having obtained the behaviour of the large deviation function of the current in the limit of small 
fluctuations, we then looked at the limit of extreme fluctuations, either positive or negative, where it turns 
out that the deformed Markov matrix can be diagonalised directly, without invoking the integrability of 
the system. In the large negative fluctuation limit, where the current goes to 0, the deformed Markov 
matrix is a perturbation of a diagonal matrix, and we find that the large deviation function of the current 
is independent of the size of the system. In the large positive fluctuation limit, where the current goes 
to infinity, the system is equivalent to an open XX spin chain, and the large deviation function of the 
current is proportional to the size of the system. We also note that these properties do not depend on 
many details of the system, and would be equally valid with added interactions between the particle or 
site-wise disorder of the jumping rates. The difference in scaling with respect to the size of the system 
suggests that a generic dynamical phase transition exists at an intermediate current for all those systems. 

Finally, we considered a less exact but more powerful approach, based on the macroscopic fluctuation 
theory, to obtain the large deviation function of any finite fluctuation of the current, using a hydro- 
dynamic description of the system in the large size limit and a non-rigorous exchange of limits. This 
allowed us to build the complete dynamical phase diagram of the current for the open ASEP. In four of 
the five phases we obtained, which are the ones corresponding to a current smaller that ^, we found large 
deviation functions independent of the size of the system, and perfect agreement with the exact results 
obtained before. In the last phase, we found no such agreement, and we surmised that the MET breaks 
down due to the necessary presence of local correlations in the system in that phase. In the phases where 
the MET does apply, we gave some evidence for a conjecture stating that the non-optimal extremisers 
of the MET action are the slow relaxation modes of the system. 

Many challenges related to the contents of this review remain to be tackled, and we will conclude by 
mentioning a few of them. 

First of all, the Q-operator method summarised in section |IV B 2| has only yielded the dominant 
eigenvalue of the deformed Markov matrix so far, but it allows in principle access to the whole spectrum. 
Obtaining it, through that method or any of the other ones mentioned in section [iVB[ would in particular 
allow to verify that the non-optimal profiles mentioned in section |VI C| are really related to the low 
excitations of the model. 

On the subject of those local minima of the MET action, it would be interesting to understand for which 
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models they might appear and whether they always correspond to excited states. For one-dimensional 
lattice gases, a finite drive in the bulk seems essential: those local minima appear for the ASEP because 
of the appearance of shocks and anti-shocks which can be combined in many ways, which is a consequence 
of the diffusion part of the Langevin equation being inversely proportional to the size of the system. It 
is natural to wonder what might be the case in other types of systems, such as higher-dimensional ones, 
for instance. 

Finally, it would certainly be interesting to understand exactly how generic the dynamical phase 
transition we discussed in section |V C| really is, and whether it has any universal features. In the case 
of the open ASEP, that transition arises because of the fact that the current which can be obtained 
from a hydrodynamic description is bounded. If one tries to impose a current higher that ^ (which 
is the upper bound), which is not forbidden in the underlying microscopic model, the hydrodynamic 
description breaks down, and another macroscopic description of the model is needed. It would stand 
to reason that such a transition would appear for any model where a hydrodynamic limit introduces 
bounds for an observable. The question of whether those transitions have any universal features is for 
the moment entirely open. 
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